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ABSTRACT 

As a consequence of the ever-shrinking sizes of nanoelectronic devices, hitherto 
neglected quantum effects, such as tunneling, are becoming important for device 
characterization The study of electron reflection and transmission probabilities at 
potential barriers is one of the important areas of active research in this field. 

Analytic solutions for the quantum-mechanical transmission coefficient through a 
potential energy profile of arbitrary shape do not exist. One conceivable method for 
finding the transmission coefficient through such a potential involves transfer matrices. 
This technique is numerically limited, unfortunately, and fails to provide adequate results 
for potentials of interest in the development of practical nanoelectronic devices. 

However, within its capabilities, the transfer matrix method is a useful reference to which 
other results may be compared. Another method, utilizing backward recurrence, has been 
proposed as a numerically stable alternative for calculating the transmission coefficient 
through such potentials. This second method has yet to be widely applied. 

This thesis investigates the capabilities and limitations of each method, with an 
emphasis on their scope of applicability. Extensive programming, in the C language, has 
been done to examine the two methods. Output from these programs has been analyzed, 
and the backward-recurrence method has been shown to have wider applicability, and to 
be faster and much more numerically stable. 
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I. INTRODUCTION 



Numerical methods are required to predict the characteristics of nanoelectronic devices, 
which rely on quantum tunneling for their operation. Continuing scale reduction of modem 
semiconductor electronics will inevitably result in devices which rely on the principles of quantum 
physics. 

In this thesis, I shall apply two different numerical methods to the problem of calculating 
the transmission coefficient. One, which relies on transfer matrices, is known to have limited 
applicability due to inherent numerical instabilities, but provides a useful benchmark in cases for 
which it is accurate. The other is a backward-recurrence algorithm, thought to be much more 
numerically stable, which potentially has wide future applicability. I have written several original 
computer programs, in the C language, which utilize both the transfer matrix and backward- 
recurrence techniques. These programs, and analysis of their output, are my contribution to this 
research. To my knowledge, the backward-recurrence technique applied herein has not been 
widely used; it proves to be powerful in solving for the transmission coefficient through difficult 
potentials. 

The values for the transmission coefficient obtained by each method are first compared to 
a known correct result: the transmission coefficient through a single, square potential energy 
barrier. In addition to this test for accuracy, the transfer matrix method contains a built-in test for 
numerical instability, which also is exploited. Examples, which reveal the numerical limitations of 
each method, are included as well. Graphical program output is analyzed to critique the 
capabilities of the two numerical methods. 
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II. QUANTUM TUNNELING 



In classical physics, the motion of a particle is governed by conservation of energy: a 
particle of total energy E is unable to move past any point beyond which the potential energy is 
greater than E. Such a point is said to comprise a potential barrier to the particle’s passage. If 
electrons in semiconductor devices truly behaved as classical particles, then the tunneling diode 
and zener diode would not exist, and aluminum household wiring could not conduct current 
through twisted-wire junctions covered by an oxide layer which, in bulk form, is an excellent 
insulator. These items all depend upon the concept of quantum tunneling. 

Quantum theory allows for the penetration of barriers of energy greater than the incident 
particle's energy. This phenomenon is a consequence of uncertainty in the position of the 
electron, according to the uncertainty principle. (Eisberg & Resnick, 1985, pp. 65-68) The 
electron's wave function, has a non-zero magnitude at all points in space, though it 

decreases rapidly beyond potential barriers. Since the probability of finding the electron in a 
specified region is proportional to the square of the magnitude of the wave function, there results 
a finite probability that the electron will be found on the other side of the barrier, in a classically 
forbidden region. 

For non-relativistic systems, motion is governed under quantum theory by the Schrodinger 
wave equation. (Eisberg & Resnick, 1985, pp. 177-209) The general three-dimensional form of 
this equation (the time-dependent Schrodinger wave equation) is 

ih — = V 2v F+F(:c,.y,r)'F. (1) 

dt 2m 
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If the potential of interest is essentially one-dimensional in form. Equation 1 simplifies to 



^ n 2 d 2 v ... ^ 

4 ¥*-SF tW <2) 

In the above equations, h=h/2tt, where h is Planck’s constant [6.626x1 O'" 1 J sec], / is the square 
root of -1, m is the mass of the particle, V(x) is the potential function, x is the spatial coordinate, 
and / is time. In the development of this equation, several key hypotheses must be made, which 
are worth reviewing. 

First, the wave equation must be consistent with de Broglie’s postulate regarding the 
wave-particle duality of matter, namely 

'4 (3) 



where p is the momentum of a particle, and A is the wavelength of the wave packet representation 
of the particle. (Thornton & Rex, 1993, pp. 178-179) Second, the wave equation must be 
consistent with the relation 



E = £-+V, (4) 

2m 

the non-relativistic expression for total energy of a particle as the sum of its kinetic and potential 
energies. 

It should be noted that, although Equation 4 holds only for particles moving at non- 
relativistic speeds, experience has shown that the Schrodinger formalism is adequate to describe 
the behavior of the systems studied in this thesis research. (Eisberg & Resnick, 1985, pp. 128- 
132) 

The form of the wave function, which is the solution to Equation 2, must next be 



4 



determined. For a free particle, one form of this solution is 

T^x,/) = cos(Ax-cy/)+/sin(£x-<y/), (5) 

where k is the wave number [27 t/X], co is the angular frequency [2nf], t is time, and /is the 
frequency. This, as required above for the case of V constant, is a traveling- wave solution. It is 
also complex-valued, which is critical to the remainder of this discussion. The physical 
significance of this complex-valued wave function, postulated by Born in 1926, is that the 
magnitude squared of the wave function is the probability density of the particle, or, stated 
another way, that 

'¥\x,t)'¥(x,t)dx (6) 

is the probability of finding the particle between x and x+dx, at time t. ¥* (x,t) here is the complex 
conjugate of the wave function. 

If the potential of interest is not a function of time, so that the system’s energy is 
conserved, then the left-hand side of Equation 2 may be associated with the characteristic 

ay 

equation ih = EEV , where E is the energy eigenvalue (a constant) of the system. The spatial 

dt 

and temporal dependencies of ¥ may then be separated, so that ¥(x,t) = iJ/( x)e' Et/h . 

The result of this separation of variables is that the time-independent Schrodinger wave 
equation is given by 

-T“T T-+r(*M*) = £>(*), (7) 

2m ax 

where V(x ) is the potential and E is the total energy. t//(x), here, represents only the spatial 
dependence (on x, in this case) of the wave function ¥(x,t). The full derivation is given by Eisberg 
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and Resnick. (Eisberg & Resnick, 1985, pp. 151-167) 

Equation 7 is an ordinary differential equation, not a partial differential equation; this is a 
major reduction in the complexity of the problem. We shall now focus on the solutions 
(eigenfunctions), y/(x), of Equation 7, which hold for a region of constant potential V: 

= Ae ,kx + Be-'*, for E>V 

YK > (g) 

y/(x) = Ce “ + De", for E <V, 
where A, B, C, and D are constants, and k is equal to /7r, with 

JlmiV - E) 

k = ±-±- (9) 

n 

Equation 9 provides the correct propagation constant only when V is constant. While 
approximating V(x) as a constant may seem limiting, it will be shown that one can obtain useful 
results in spite of this approximation. Note the presence of the increasing exponential e K in 
Equation 8. This observation will later prove critical. 

The boundary conditions for the spatial wave function are that y/(x), and its first derivative 
with respect to x, must be finite, single-valued, and continuous. Eisberg and Resnick’s description 
of the reasons for these restrictions is excellent. (Eisberg & Resnick, 1985, pp. 155-157) We will 
utilize these boundary conditions in what follows. 

As an example, let us apply these solutions to the case of a rectangular potential barrier. 
Let this barrier potential be V(x)=0 for x<-a and x>a, V(x)=V 0 for -a<x<a. Classically, if the total 
energy of a particle is such that E> V Q , then the particle always passes the barrier; similarly, if 
E<V 0 , then the particle never passes the barrier. Since the spatial wave function y/(x) has the 
form of Equation 8, for the quantum mechanical case, the behavior of the particle is different. 
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First, if E> V Q , there will be a nonzero probability of reflection, at both edges of the barrier. 

Stated differently, there is less than a 100% chance that it will traverse the barrier. More 
importantly, however, if E<~ V 0 , there will be a nonzero probability of transmission through the 
barrier, or tunneling, also known as barrier penetration. 

The form of y/(x) in the region x< -a is Ae lkx +Be' lkx , corresponding to a pair of incoming 
and outgoing traveling waves. The form of the general solution in the region x>a is similar, 
namely Fe ,kx +Ge ,kx . Between -a and a , however, the form is Ce^+De**. Boundary conditions 
require that, at x--a, Ae ,ka +Be‘ ka =Ce Ka +De Ka so that i//(x) is continuous. In order for the first 
derivative d y/(x)/dx to be continuous, A e ,ka -Be‘ ka = ( i K'k)[Ce Ka -De Ka ] . It is shown, in Appendix A, 
that the amplitudes A and B are related to the amplitudes F and G by a matrix involving k, k, and 
the barrier width a. This matrix is called the overall transfer matrix, M. /rand k , of course, 
depend directly on the particle mass m, the total particle energy £, and the barrier height V 0 , as 
shown previously. In short, there are four physical parameters which determine M: particle mass, 
energy, barrier height, and barrier width. The form of the overall transfer matrix for the case of a 
simple rectangular potential is given in section A of Appendix A. The overall transfer matrix has 
some interesting symmetry properties; these are developed in Appendix A. The determinant of 
the overall transfer matrix in also discussed in Appendix A. 
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III. THE NUMERICAL METHODS 



A. OVERVIEW 

Quantum-mechanical transmission through a potential barrier is of interest. Current 
microprocessors have a device scale of 0.3 pm. The room-temperature de Broglie wavelength of 
a conduction electron in silicon is approximately 0.025 pm. Future device scale reductions, by 
merely a factor of ten, will therefore introduce significant quantum behavior. If the current pace 
of microprocessor scale reduction continues, knowledge of the quantum-mechanical transmission 
coefficient will be required in the design of electronic devices within the next 1 5 years. 
(Semiconductor Industry Association, 1994) The physical meaning of the transmission coefficient 
is described in full by Eisberg and Resnick. (Eisberg & Resnick, 1985, pp. 193-198) Essentially, 
the transmission coefficient is a measure of the probability of barrier penetration by a quantum 
particle, such as an electron. 

In the preceding chapter, we discussed the wave function for a free particle. A conduction 
electron in a semiconductor, however, is not a free particle: it experiences the spatially periodic 
potential energy environment of the crystalline lattice. Here we invoke the effective-mass 
approximation, which, for electron energies near the band energy extrema, permits us to treat the 
motion of the electron in the solid as if it were a free particle having an effective mass, m *. The 
value of the effective mass is a material-specific parameter. For most devices of technological 
interest, the electron energies remain near the band extrema and the effective-mass approximation 
is valid. A complete description of the effective mass is given by Eisberg and Resnick. (Eisberg 
& Resnick, 1985, pp. 460-464) 
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The transmission coefficient is defined as 



T = 



F 
A ’ 



( 10 ) 



where the coefficients F and A are the magnitudes of the transmitted wave function and the 
incident wave function, respectively. Let us consider a barrier which extends from x=xl to x=xr. 
The potential outside the barrier is constant. The incident and reflected wave functions have the 
form 



V', naden M = Ae‘ k ' x 

V reflect) =Be" k '\ 

for x<Xi, and the transmitted wave function is of the form 

W transmitted (^) = , 

for x>Xr, where k t is the wave number in the incident medium (before the barrier) and k 2 is the 
wave number in the transmitted medium (after the barrier). The wave number is a complex- 
valued quantity defined by k=iK, and at is given by equation 9. Substituting the effective mass, 
m *, for the mass m , results in the equation 



(ID 



( 12 ) 



1 2 tn{E-V) 

\ lc~ 



(13) 



where V is the constant value of the potential energy. 

The two numerical methods for finding the transmission coefficient described herein both 
use solutions based upon piecewise-constant potentials. The first numerical method is the transfer 
matrix method, which starts by breaking up the actual potential into an approximate staircase 
potential, which is piecewise constant. A two-by-two transfer matrix (N matrix) is then calculated 



10 



at each interface, and an overall transfer matrix (M matrix) for the entire potential is found by 
matrix multiplication of the individual interfaces' matrices. The transmission coefficient for the 
entire potential can then be found. A more detailed description of this procedure will be given 
later. 

The second method is a backward-recurrence algorithm, suggested by Luscombe 
(Luscombe, 1992, pp. 1-20), which greatly increases numerical stability by eliminating the error 
associated with the forward propagation of a solution involving both decaying and growing 
exponential components. A complete description of this technique also appears later. 

This thesis emphasizes the applications and numerical accuracy of these two techniques. 

It will be shown that the transfer matrix technique has a convenient built-in test for numerical 
instability. The mechanism of the numerical instability will be discussed in detail, as will the 
conditions which bring it about 

The backward-recurrence method has excellent numerical stability, and also runs faster. It 
does not have a comparable built-in method to test for accuracy, however. The results of 
applying the method to complex potentials will be discussed. 

There exists a closed-form analytical solution to the simple tunneling problem, for the case 
of one square barrier in a region of otherwise constant potential. (Singh, 1997, pp. 131-135) To 
test the accuracy of both numerical methods, each technique’s accuracy in giving a numerical 
solution to this known problem was checked. 

All computer programs used in this thesis were written in ANSI-standard C. This 
language allows the specification of either single- or double-precision floating point numbers, and 
thus provides the ability to investigate numerical and algorithmic stability at two different 
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precision levels. The MATLAB programming language (which has built-in matrix functionality) 
was not used in this work, primarily because it lacks this capability, but also due to its slower 
operation in situations requiring loops, which prove to be unavoidable in finding the transmission 
coefficient at multiple values of incident particle energy. 

B. TRANSFER MATRIX METHOD 

The interface transfer matrices were the subject of a particularly concise treatment by 
Singh, which 1 paraphrase here: (Singh, 1997, pp. 148-149) 

The transfer matrix method assumes a piecewise-constant potential approximation to the 
potential of interest. Therefore, the wave function has the form 

y(x) = Ae ,kx +Be- ,la (14) 

in each region of constant potential, as shown earlier The incident wave is assumed to be 
travelling in the positive x-direction Each region will have its own set of coefficients A and B and 
wave numbers k. Assuming that the electron's effective mass is constant in all regions, and again 
applying the boundary conditions at the interface between regions, here situated at x=x 0 , 

HO = KO 

dy/ _dy/ 05) 
dx x ° dx x °* 

we see that the coefficients A and B can be related to each other by the matrix equation 



Ul 




fAl 




= N 


L 


A. 




A2. 



where A / and Bj are the coefficients of the wave function for x x 0 , and A? and B 2 are the 
corresponding coefficients for x- x 0 . The matrix N is called the interface transfer matrix. This 
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two-by-two matrix can be shown, in the general case of a potential with interfaces at points Xj 
separating regions of potential F, +/ , to have components 



AW) 

N,(2,2) 



]_ 

2 

~2 

s 



1 +- 



iC; 



»+i 



j _ k, + \ )(*,♦!+*<) 



1--^ 






V < J 



1 + 






(17) 



The matrix product of all of the N matrices yields the total transfer matrix for the potential, M. 
We use the total transfer matrix to calculate the overall transmission coefficient, by the relation 



|A/(U)|‘ 



(18) 



This method is robust for simple potentials, but runs into difficulties if the barrier height, barrier 
width, or number of interfaces is too great. This will be discussed in greater detail later. The 
program shows perfect agreement with the analytic form of the transmission coefficient for one 
square barrier. 

C. BACKWARD-RECURRENCE METHOD 

The backward-recurrence method used in this thesis is that described by Luscombe. 

(Luscombe, 1992, pp. 1-20) As he states, 

relevant device variables must ... be specified within close tolerances, and clearly 
the most efficient means of developing a realistic [microprocessor] design is 
through computer modeling based on fundamental physical laws, e.g. the Poisson 
and Schrodinger equations. By providing the tools to explore, in detail, and 
before fabrication, the effects of varying ... [device] parameters ..., modeling 
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optimizes the development cycle and provides a cost-effective method for 
validating and refining device concepts. (Luscombe, 1992, p.2) 

Before discussing the backward-recurrence method in detail, we note the following 
observations. The Schrodinger equation is a second-order differential equation, and, as such, has 
two linearly independent solutions. In the classically forbidden potential regions, the physical 
solution for the wave function ij/(x) is strictly decreasing as a function of distance. (Since E<V in 
the classically forbidden regions, the second derivative of iff is strictly positive.) Now, as is simple 
to check, for example by examining the Wronskian relation associated with the Schrodinger 
equation, if one solution is monotonically decreasing, the other, linearly independent solution is 
monotonically increasing. This second, linearly independent solution is therefore not physical. In 
numerical methods that propagate, through the use of iteration, both the growing and decaying 
solutions on equal footing (as in the transfer matrix method), the slightest round-off error will 
trigger the growth of the unwanted, growing solution in the classically forbidden potential 
regions. By employing backward iteration, however, the physical solution becomes dominant, 
increasing resolution in the direction of backward iteration. 

The transfer matrix method, while accurate in limited applications, is numerically unstable, 
as stated above. This instability is implicit in the N matrices (interface matrices) described earlier, 
(discussions with James Luscombe, May-June 1997) 

The backward-recurrence method works as follows. We shall assume the electron's 
effective mass to be constant throughout, and equal to 

m =0.061 m e \ (19) 
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we therefore need not change it from region to region {m e , here, is the electron mass). The only 
thing that changes spatially is the potential, V. As with any numerical solution to a differential 
equation, we employ the finite difference approximation. We sample the wave function at discrete 
points x„=nA, where A is the step size, and is chosen to be sufficiently small to accurately resolve 
the potential. We begin with the incremental Taylor series expansions 



t//(x + A) = y(x) + A y' (x) + — y" (x) + {higher - order terms) 

A 2 

y[x - A) = y/(x) - A y'{x) + — y" (x) + {higher - order terms). 



( 20 ) 



An expression for y”. 



y . (x)= H* + A) + lK*-A)-2^ (2() 



can then be derived. The algorithm progresses from right to left, through the perturbation 
described below, in spatial steps of size A. The solution at the / 1 th step may be written as 



Vn = 



( 22 ) 



The time-independent Schrodinger equation. 



ft 2 



- — r y"+V y = Ey, (23) 

2m 



is equivalent to a three-term recurrence relation, 

^ n+l +^ n -,+^ n = 0 ^ ( 24 ) 



where h„ is given by 



K=mf\{E-V n )-l. (25) 
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Finally, using the definition 




the above three-term recurrence relation can be rewritten as 



which is easily recognizable as a two-term backward-recurrence relation It is this relation which 
is the heart of the backward-recurrence method. 

The potential being analyzed must consist of an irregularity of finite spatial dimensions, 
beyond which the potential is constant. The irregularity in the potential may be of arbitrary height 
(eV), curvature, and width (nm). There must also exist two regions of constant potential, which 
we designate as right-hand and left-hand regions. Each has an associated wave number, k R and k L . 
The constant potential energies in these two regions need not be the same, necessarily, but for 
simplicity they are assumed to be identical in this thesis. 

Defining r Q to be the value of r at the edge of the left-hand region, it is easily seen that 



Similar logic can be applied to the right-hand region of constant potential, by writing r N+ i as 



t 




( 28 ) 



where iff„=ilr(x)=e kLx + e' ,kLx and x=nA, and that the reflection coefficient R is given by 
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- II. 

'N + \ ~ 



i*j,A(Af+l) 



Vtt Te lk « A(N) 
where y/ ys(x) = Te ,k * x . 



= e ,k * A = r, 



N ■> 



(30) 



We know r N+t if we know k R and the step size, A. The expression for k R is 



k ’=l¥) E - y *>- (3i) 



where V R is the constant potential in the right-hand region. An analogous expression exists for k L , 



k, = 



2m 



(E-FJ. (32) 



We have derived one expression for R, Equation 29. T is related to R by the expression 




sin k L A 
sin k R A 




(33) 



where T and R are the overall transmission and reflection coefficients, respectively. (Luscombe, 
1992, p. 7) Equation 33 may be rewritten as 




^sin&^A^ 


(4 sin£ L A)lm(r 0 ) ^ 


^sin k R A J 


[l-2Re( V “‘ 1 )+|r.| 2 J 



(34) 



which the program uses to calculate T once it has found the value of r 0 . 



D. ANALYTIC SOLUTION 

The analytic solution to the problem of one square barrier, used to test the transfer matrix 
and backward-recurrence methods, is simple to derive. (Eisberg & Resnick, 1985, pp. 199-201) 
(Singh, 1997, pp. 131-134) Consider a square barrier of height V 0 and width a, extending from 
x=0 to x=a. The potential everywhere outside the barrier is zero. Define the region to the left 
and right of the barrier to have wave number kj, and the barrier to have associated wave number 
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kjj. These wave numbers are given by 



k, = 
k u = 




(£- 0 ) 



Ie-v.Y 



(35) 



As stated earlier, the solution to the Schrodinger wave equation in a region of constant potential 

has the form y/(x) = Ae lkx + Be ~ ,kx . Call the region in which x- 0 region I, the region in which 

0<x<a region II, and the region in which x a region III. The wave functions in these regions are 

y/, = Ae‘ k ' x +Be- ,k ' x 
y/ n =Ce‘ k " x +De-' kuX (36) 

Vm =te' k ' X , 



assuming that there is only a right-moving wave in region III. Applying the boundary conditions 
on the wave functions given by Equations 15, atx=0 and x=a , and performing algebraic 
simplification, one arrives at the expression 



F 


2 


4k, k^"-^ 


A 




(k, +k, l y ~(k, -*„)V'** 



(37) 



which is the expression plotted by the programs to check the accuracy of the two numerical 
methods. 
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IV. THE PROGRAMS 



A. TRANSFER MATRIX PROGRAM 

A copy of the transfer matrix program, method], c, is attached as Appendix C, Section A 
The program is written in ANSI C. 

Since ANSI C contains no built-in complex numbers capability, a global variable structure 
called complex is defined, consisting of two double-precision floating point numbers representing 
the real and imaginary parts of the complex number. Also, since the potential, V, and wave 
number, k, as functions of x, are required in multiple functions, they are defined globally. All 
other variables have local scope. 

The functions addc, subc, mule, and dive have been written to perform the four basic 
mathematical functions on complex arguments. They all take two arguments of type complex, 
and return type complex. The function expc computes where x is a complex argument, and 
returns type complex. The function absc, given a complex argument, returns the (real) magnitude 
of the complex number, as a double-precision floating-point number. The function 
determinantcomplex takes four arguments of type complex, which represent the elements of a 
two-by-two overall transfer matrix, and returns the (complex) determinant of the matrix. 

N1 1 complex, N12complex, N21 complex, and N22complex calculate the values of the 
elements of the interface transfer matrix, using Equation 17. They take four arguments of type 
complex, representing the wave vector left and right of the interface, their ratio, and the value of x 
at the interface. 
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Finally, the function makeV initializes the potential energy array, V. In this program, the 
potential energy profile is a series of square barriers of height Vo, width BARWIDTH , and 
separated by BARWIDTH. BARRIERS represents the number of barriers present. ERIGHT, 
ELEFT, and EPTS represent the highest and lowest values of incident particle energy, and the 
number of energy values used, respectively. Vo, BARWIDTH, BARRIERS, ERIGHT, ELEFT, 
and EPTS are all predefined constants, listed as #define statements. 

Three other physical constants are used. One is the effective mass, defined as a ratio, 

* 

ftl 

— , equal to 0.067, and represented by EFFMASS. This is an appropriate value of the effective 
m e 

h 2 2 

mass for GaAs semiconductors. The value of , in units of eVnm , is defined by H201ER2M 

2m e 

as 0.0381, and is constant for all materials. Finally, the program calculates a constant, C, which is 

, 2 m _ L . , , J EFFMASS 

equal to — — . C, obviously, may be expressed as . 

h 1 H20VER2M 

Function main generates the program's output, and controls the execution of the transfer 
matrix method algorithm. The printf("%. 12f\t%. 12f\n",E,T); statement lists the values of incident 
particle energy and transmission coefficient in two columns, to the screen. To generate the graphs 
found herein, the program's output has been saved in a file using indirection, and plotted using the 
commercially available program Spyglass Plot. 

The basic execution of main is as follows. A value of incident particle energy is assigned 
by the for(ec=0.0;ec<epts;ec+=1.0) statement. This loop goes through all desired values of 
energy, from ERIGHT down to ELEFT, dividing the energies into EPTS steps. As EPTS is an 
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integer, the variable epts is introduced to promote it to a double-precision floating-point number 
in order to be type compatible with the variable ec. This also prevents errors, which might be 
caused by inadvertent integer division. 

Once a value of E has been defined by the E=ERIGHT-(ERlGHT-ELEFT)*(ec/epts); 
statement, the variables Mil, M12 , M21 , and M22 , all of type complex, are initialized. These 
variables are the matrix elements of the overall transfer matrix. Since the overall transfer matrix 
will be calculated by matrix multiplication of all intervening interface transfer matrices, it is 
initialized as a two-by-two identity matrix by the 

oldM 1 1 . imag=oldM 1 2 . imaged dM2 1 . imag=oldM22.imag=0 . 0; 

oldMl 1 . real=oldM22 . real= 1 .0; 

oldM12.real=oldM21. real-0.0; 

statements. 

The next step is to initialize the wave vector array, k, at the current value of E. The 
for(xc=0;xc<XPTS;xc++) loop handles this. At each point, the wave vector will be either purely 
real or purely imaginary, and is given by Equation 13. 

The for(xc=l;xc<XPTS;xc++) loop, which follows, goes through the potential profile and 
identifies interfaces; that is, it finds values of x at which k changes. When an interface is found, 
the if((k[xc].real!=k[xc-l].real)||(k[xc].imag!=k[xc-l].imag)) statement executes, causing the 
interface transfer matrix to be calculated and matrix multiplied by the "running total" overall 
transfer matrix, given by elements oldMl 1, oldM12, oldM21 , and oldM22, all of type complex. 
Upon completion of this if statement, the program checks for numerical instability. 

The Mdet=determinantcomplex(newM 1 l,newM12,newM21,newM22); statement checks 
the determinant of the calculated overall transfer matrix. The subsequent 
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if((Mdet.real>1.0001)||(Mdet.real<0.9999)) statement checks to see if the determinant is one. If it 
is not, within four-digit precision, then the program halts execution and prints the listed error 
message. 

The transmission coefficient, T, is then found using Equation 1 8, by the 
T=sqrt(1.0/(absc(newMl l)*absc(newMl 1))); statement. It is then printed out in two columns, as 
mentioned before. Spyglass Plot then produces graphs from the output. 

B. BACKWARD-RECURRENCE PROGRAM 

A copy of method2.c, which employs the backward-recurrence method, also is included 
in Appendix C, Section B. The functions addc, subc, mule, dive, expc, and absc are re-used here, 
to perform the basic operations on complex numbers. The predefined constants are also the same 
as those used in method l.c, with the exception that VL, VM , and VR have been added to allow the 
specification of different potential in the left-hand and right-hand constant-potential regions, and 
between the potential barriers. The program operates as follows: 

Variables r, ro , rnplusl , and bnfXPTS ] , of type complex, hold the values required in the 
backward-recurrence relation given by Equation 27. Function makeV, as before, creates a 
potential profile consisting of square barriers. The barriers have height Vo, and separation and 
width given by BARWIDTH. Other versions of this program have been written which calculate 
the transmission coefficient through different potentials. Only makeV need be changed to 
accomplish this. 

The for(ec=0.0;ec<epts;ec+=1.0) loop counts through the values of incident particle 
energy, from BRIGHT down to ELEFT, using the E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); 
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statement. At each value of incident particle energy, the for(xc=0;xc<XPTS;xc-H-) loop initializes 

the bn array, bn is given by Equation 25, which is 

bn[xc] . real=C * delta* del ta*(E-V[xc] )-2 . 0; 
bn[xc].imag=0.0; 

in the C language. 

kl and kr represent the wave numbers in the left-hand and right-hand regions of constant 

potential, respectively, and are given by Equations 31 and 32 and by 

kl=sqrt(C*(E-VL)); 
ki=sqrt(C*(E-VR)); . 

Equation 30 gives rv+/, the value of r in the right-hand region of constant potential. The variable 

rnplusl represents r^+i in the program, and it is initialized by the 

mplusl ,real=cos(kr* delta); 
rnplusl .imag=sin(kr*delta); 



statements. 

The for(xc=(XPTS-BARPTS);xc>=BARPTS;xc— ) loop, while deceptively short, actually 

performs all of the backward-recurrence steps, using Equation 27, and updating the value of 

rnplusl before each iteration. Finally, having computed ro , the program calculates the value of 

the transmission coefficient, T, using Equation 34. This is done by the 

temp.real=0.0; 

temp.imag=kl*delta; 

temp2=mulc(ro,expc(temp)); 

T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(l .0- 
2.0*temp2.real+absc(ro)*absc(ro))); 
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statements. The values of E and T are then printed out, in columns, as before, and the results 



plotted using Spyglass Plot. 
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V. PERFORMANCE AND NUMERICAL STABILITY 



A. TRANSFER MATRIX METHOD 

The transfer matrix method, as mentioned, contains a testable quantity which allows one 
to determine when its output has begun to be unreliable. It is shown in Section B of Appendix A 
that the determinant of the overall transfer matrix, M, for a system must equal one. This has been 
used to find the point at which the program's output begins to deteriorate. 

By testing the determinant in this way, one also may find the specific value of incident 
particle energy at which the transfer matrix method breaks down. Extensive results of these tests 
are attached as Appendix B, for a simple test potential composed of a series of square barriers, 
whose number, height, and width/separation may be changed until the program detects numerical 
instability. From these test cases, it is clear that the transfer matrix method can be numerically 
unstable in many applications. Examples of this instability will be given shortly. 

In Appendix B, Section A is the output of the transfer matrix method program, 
method l.c, when single-precision floating-point numbers are used. These numbers have the 
equivalent of six decimal places of precision. Consider the case in Appendix B, Section A, Part 1 
(B. A. 1 .) , of one five nanometer (nm) barrier which is 0.2 electron volts (eV) tall. In this case, for 
single-precision numbers, the method breaks down at an incident particle energy of 0.010702 eV. 

method l.c analyzes particle energies beginning with the largest value, so for any particle energy 
less than this value, the transfer matrix method returns inaccurate results. Compare this result to 
that in B.B. 1., which is the same physical case, analyzed using double-precision floating-point 
numbers, which carry 12 digits of precision. No numerical instability occurs with the double- 
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precision numbers; this is not surprising, since more precision is available in the calculations used 
in the program. There is no number of digits of precision, however, which is sufficient to entirely 
prevent numerical instability of the transfer matrix method. Note in B.B.I., for instance, that 
inaccuracy occurs when the height of a single, 5 nm barrier is increased to 5 eV. 

In B.A.2. and B.B.2., the barrier width has been varied. Note that for both single- and 
double-precision floating-point numbers, there exists a value of barrier width which causes 
numerical instability of the transfer matrix method. Similarly, in B.A.3. and B.B.3., this is 
observed when the number of square potential energy barriers is increased. In each case, the 
value of incident particle energy at which instability occurs is listed. Observe that, as the 
parameters of the calculation are pushed beyond the point at which instability just starts, the 
minimum incident particle energy rises. This effect can be seen for either single- or double- 
precision floating-point numbers, throughout Appendix B. 

The question arises: how severe is this effect in the transfer matrix method? At this point, 
method l.c is applied to a situation with a known, analytic solution, in order to find out the answer 
to this question. The following paragraphs will refer to various figures, all of which have been 
collated in Appendix D. 

The test case chosen is that of one square potential energy barrier, of width 5.0 nm and 
height 0.23 eV. In B.A.2., it has been shown that the transfer matrix method breaks down at 
incident particle energy 0.013768 eV, for single-precision numbers. It does not fail if double- 
precision numbers are used. Figure 1 shows the output of the program using double-precision 
numbers This figure is a familiar sight in introductory quantum physics texts. (Eisberg & 
Resnick, 1985, p. 202) (Singh, 1997, p. 135) There is a positive transmission coefficient for 



26 



incident particle energies less than the barrier height, indicating that quantum tunneling will occur. 
There is also a region of non-unity transmission coefficient at values above the barrier height, 
indicating the presence of quantum scattering. 

In Figure 2, the transfer matrix method solution minus the value of the transmission 
coefficient calculated by the analytic solution for one square barrier has been plotted. Note that 
the transfer matrix solution is exact. The same figure resulted from running the program with 
single-precision floating-point numbers. It is apparent that, when the determinant of the overall 
transfer matrix is 1.0001 instead of 1.0000, significant degradation of the numerical solution is not 
observed. 

This is an interesting dilemma, since analytic solutions for the transmission coefficient are 
so rare in practical quantum tunneling problems. When is the output of the transfer matrix 
method useful? What are the limitations of the method? Only qualitative answers are available to 
the first question; the method is useful when it is numerically stable. Instability may be detected 
by means of the determinant of the overall transfer matrix. The limitations of the transfer matrix 
method, however, are tested in this thesis by five test potentials; a single square barrier, sequences 
of two, three, and five parabolic barriers, and a resonant-tunneling diode (RTD) potential. 

The results, in the case of the RTD potential, clearly demonstrate the inadequacy of the 
transfer matrix method in dealing with elaborate potential profiles. Appendix C, Section E 
describes mlpara.c, which tests the transfer matrix method against potential energy profiles 
composed of series of parabolic barriers, rather than square barriers. An example of a potential 
used in this program is given in Figure 3. In this figure, one can clearly see the difficulty posed by 
potentials consisting of smooth curves; the potential must be carved into sub-bins, each of which 
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is piecewise-constant. The partitioning must be such that the character of the potential in 
preserved, which, for steep slopes in the potential, may require small step sizes. Recall that the 
number of bins (barriers) tends to drive the method into instability, as shown previously. 

Figures 4, 5, and 6 are the output of mlpara.c for the cases of 2, 3, and 5 parabolic 
barriers. In each of these cases, the “safety check” of the determinant of the overall transfer 
matrix has been performed for each value of energy, and program execution halted when the 
determinant differs from 1.0 by more than 0.0001. In the case shown in Figure 4, instability did 
not occur. In Figure 5, it occurred after 952 out of 1000 energy points had been analyzed. In 
Figure 6, it happened after only 908. Each graph, however, shows only the values of T which 
were calculated before instability occurred. These values are therefore known to be accurate. 

It has been established that the output of the transfer matrix method is correct, even when 
the determinant of the overall transfer matrix differs slightly from one. When the discrepancy in 
the determinant is less than 0.0001, therefore, the prpgram’s output is accurate. These figures 
will be used as a basis for evaluating the performance of the backward-recurrence method, applied 
to the same potentials, in Section B of this chapter. 

The results seen in Figure 4 seem reasonable. There is a peak transmission coefficient of 
1.0 when the energy equals the height of the barriers, which agrees with observations based on 
systems of square barriers. Scattering occurs for incident energies above the barrier energy, as 
expected. Since there are two barriers, one expects a resonant-tunneling energy to exist, and 
there is, in fact, a peak in the T versus E curve for incident energies below the barrier energy. The 
performance of the transfer matrix method in this case is acceptable. 
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Figure 5, for the case of three parabolic barriers, demonstrates that the technique is fairly 
accurate, as well. There is a scattering region above the barrier energy, as expected, and there is 
also a resonant-tunneling energy. In addition, we see the expected twinning of transmission 
resonances due to the presence of two energy wells between the three barriers. However, as noted 
above, the method breaks down at low energies. 

In Figure 6, similar performance occurs. This figure is plotted with a logarithmic vertical 
axis, to better show the range of values. Again, twinning of transmission resonances appears, 
with four resonances this time. This is due to the four energy wells that exist between the five 
energy barriers. Performance degrades at low energies once more, and instability sets in at lower 
energies this time, as anticipated. 

Finally, the transfer matrix method is applied to a truly difficult (and physically relevant) 
potential: the resonant tunneling diode (RTD) potential. This potential is shown in Figure 7, for 
an AlAs/InGaAs/InAs RTD. (Luscombe, 1992, p. 4) Note the steeply-sloping sides of the central 
well, in combination with the smooth curve on which the central well rests. These features make 
it extremely difficult to approximate this potential by a piecewise-constant model. 

Figure 8 shows the results obtained from the transfer matrix method in this case. The 
program’s output, even using double-precision numbers, is so inaccurate that the transmission 
coefficient is not calculated between 0.75 and 0.87 eV, and the method breaks down completely 
below about 0.55 eV. Only the energies of the two peaks are correctly predicted, and the method 
completely misses a third, low-energy peak, as described in Section B of this chapter. 

In defense of the method, it must be said that for those potentials which are simple enough 
for it to handle, it performs perfectly. However, the transfer matrix method does have the 
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observable disadvantage of considerably longer execution time, for the same potential, when 
compared to the backward-recurrence method. The observed time difference has been as much as 
several minutes, on a Sun SparcStation 5 workstation, depending on the specific problem. 

The method also proves to be able to operate accurately even when the determinant of the 
overall transfer matrix differs slightly from one, but it rapidly becomes unstable for determinant 
values which differ from one by more than 0.001 . Such differences may easily occur if the 
number of interfaces is greater than ten. Difficulties arise when the program is required to 
continue to handle many more, lower, incident particle energies, beyond that energy which first 
caused the determinant not to equal one. The incident energies are analyzed in decreasing order, 
so that the method first encounters the cases least likely to cause numerical instability. Cases do 
exist, however, for which the transfer matrix method is grossly incapable of calculating the 
transmission coefficient, and it fails completely. As these cases are of significant interest in the 
development of nanoelectronic devices (like the RTD), the transfer matrix method proves not to 
be powerful enough for these applications. 



B. BACKWARD-RECURRENCE METHOD 

The backward-recurrence method has proven to be considerably faster, easier to program, 
and less prone to numerical instability than the transfer matrix method. In fact, it has not been 
demonstrated to fail numerically at all, though the fine detail of the T versus E curves of some 
potentials can also be a factor, as will be shown below, on page 32. 
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Since it lacks a built-in test for numerical accuracy, the best tests for the backward- 
recurrence method are comparisons with known transmission coefficient versus energy curves. 

For example, again consider the case of one square barrier. 

In Figure 9, the output of the backward-recurrence program has been plotted on the same 
axis as that of the analytic expression for T for one square barrier. Note the close agreement 
between the numerical and analytic solutions. In Figure 10, the difference between the analytic 
result and the backward-recurrence result is shown. The agreement is perfect, to the limits of 
precision of the calculation. 

The only other reference which one can compare to the results of the backward-recurrence 
method are known correct results from the transfer matrix method. (This makes clear why the 
transfer matrix method has been considered: it can tell the user when it has generated correct 
results.) The T versus E profiles of Figures 4 through 6, generated using the transfer matrix 
technique, can be used for comparison. 

The case depicted in Figure 1 1 is the same as that of Figure 4, except calculated with the 
backward-recurrence technique. Comparison of the two graphs suggests that the agreement of 
these two figures is perfect. From this analysis, it is clear that the backward-recurrence results and 
the transfer matrix results agree to within the calculations’ precision. 

Similar comparisons have been done between the data of Figure 12 and Figure 5, and 
between that of Figure 13 and Figure 6, for the cases of three and five parabolic barriers. For 
both of these potentials, the backward-recurrence method exactly reproduces the known-correct 
results obtained from the transfer matrix method. The backward-recurrence method executes 
these calculations in approximately one-tenth the time required by the other algorithm, as well. 
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Plot resolution impacts the quality of the program output, even when the individual data 
pairs (£,7) are correct. The rate of change of T with changes in E is so sharp in the neighborhood 
of the maxima, that with 1000 energy points, the peak itself is stepped over. Figure 13 contains 
this error, which is unrelated to the accuracy of the individual values of T. We see in Figure 13 
that the magnitudes of T at the four maxima located at about 0.06 eV are not all 1 .0. They should 
be, as Figure 14, which uses 10,000 energy points instead of 1000, demonstrates. 

As an example of this phenomenon, compare the leftmost maximum near 0.06 eV in 
Figure 13, with the same maximum in Figure 14. In Figure 13, the frequency at which energy 
values are sampled is insufficient, and this T maximum appears to be less than 1.0. In Figure 14, 
the energy values are sampled frequently enough to show the true maximum of 1.0. Required 
graphical resolution for plots of this type may easily exceed the specified energy sample rate. This 
highlights the importance of correct sample settings to the proper use of these numerical 
techniques. Fortunately, for nanoelectronic devices, the number of layers, and thus the number of 
material interfaces, will be finite; this will alleviate some of this problem, since it will limit the 
slope of, and the number of maxima in, the T versus E profile of the device. 

Using backward recurrence, the transmission coefficient for the RTD potential profile was 
easily calculated. These results have been included as Figure 15. This figure was plotted with 
10,000 energy points after noting that the maximum at about 0.6 eV had a value of T which was 
less than 1.0, when plotted using 1000 energy points. Note the full coverage, lack of numerical 
instability at low energies, and three energy peaks. This result is in agreement with Luscombe’s 
original paper, taking into consideration differences in the effective masses of the media. 
(Luscombe, 1992, p. 4) The backward-recurrence method proves itself to be quite capable of 
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attacking the RTD potential, and it no doubt is capable of analyzing the potential profiles of other 
nanoelectronic devices. This method has great promise, and wide future applicability. 
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VI. CONCLUSION 



The backward-recurrence method has much greater stability and much faster speed of 
execution, as well as a much wider scope of applicability than does the transfer matrix approach. 
Both methods reproduce, with no observable error, the transmission coefficient (7) versus 
incident particle energy ( E) curve for the classic square potential barrier, as seen in Figures 1 and 
2. 

That the determinant of the transfer matrix equals one, means that the results produced by 
the transfer matrix method are credible. This fact has been used to calculate T versus E curves for 
other potentials, such as the parabolic barrier potentials seen in Figures 4 through 6. These T 
versus E curves found with the transfer matrix method have been compared with the output from 
the backward-recurrence method, when applied to the same potential. Regardless of the 
potential’s shape, it has been shown that the backward-recurrence method produces the same 
output as the transfer matrix method, when the potential does not cause the transfer matrix 
method to become instable. 

The backward-recurrence method shows great promise in the numerical solutions of T 
versus E curves for potential energy profiles encountered in real devices, like the resonant 
tunneling diode potential shown in Figure 7. In contrast, the transfer matrix method proves itself 
to be incapable of any reasonable accuracy in this case. Other practical potential profiles for 
nanoelectronic devices also are probably within the capabilities of this algorithm. 

The backward-recurrence method for calculating transmission coefficients has been shown 
to be worth further development. The transfer matrix method, while limited in applicability, can 
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act as a worthwhile benchmark, against which future versions of the backward-recurrence code 



may be tested. 
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APPENDIX A. 



MATHEMATICAL PROPERTIES OF TRANSFER MATRICES 

A. SYMMETRY OF THE OVERALL TRANSFER MATRIX, M 
1. The Example of One Square Barrier 

It can be shown that, for a single square barrier of height V a and width 2a, the 



the magnitude of the purely imaginary wave nurrtber inside the barrier ( k bamer = itc ). 
Also in this expression for M are two aggregate constants, rj and s , which are given by 



Since, in this construction, k and k are purely real, e and 77 are real as well. 
(Merzbacher, 1961, pp. 91-92) 

2. The General Case 

The symmetry property of the overall transfer matrix, M, is obvious from its form 
given above. M has the form 



overall transfer matrix, M = 



1,1) M( 1,2) 

M(2,l) M( 2,2) 



is given by 




where k = 
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M = 



M(\,\) A/(l,2) a, +b x i a 2 +b 2 i 

M( 2,1) M{2,2) a 3 +b 3 i a 4 +b 4 i_' 



where 02=03=0, ai=a 4 , bj= -b 4 , and bj= -b 3 . In other words, the (1,2) and (2,1) elements 
of M are complex conjugates of each other, and are both purely imaginary. The (1,1) and 
(2,2) elements are also complex conjugates of each other, but each has both real and 
imaginary components. This symmetry is not limited to the case of one square barrier; 
rather, the one square barrier case has been provided as an illustrative example. 

B. DETERMINANT OF THE OVERALL TRANSFER MATRIX, M 

The determinant of the overall transfer matrix, M, is identically one. This 
condition is true for all potentials, not just for square barriers. This property can be 
verified in the M given for the square barrier example, as follows: 



C. SYMMETRY AND PROGRESSION OF INTERFACE TRANSFER 
MATRICES, N 

The interface transfer matrices have different symmetry than does the overall 
transfer matrix, M. In fact, when the numerical version of the transfer matrix solution 
begins to degrade, it is these properties of the interface transfer matrices which are 
involved. The following is output from the transfer matrix program for a condition 
known to cause numerical instability, taken from Appendix B, section B.3. The output is 




cosh2KW+ — sinh2Kior cosh 2 m sinh2*r<7 
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for a sequence of square barriers whose up-steps occur at jc= 5, 1 5,25 . . . and whose down 
steps happen atx=10,20,30.... 

BARRIERS = 10 

BARWIDTH = 5.000000 nm 

Vo = 0.230000 eV 
E=0. 196624 eV 



When x=5. 000000, N is 
-1.3361 13e-01 +-8.990499e-02 j 
-1.3361 13e-01 +8.990499e-02 j 


-1. 78344 le+00 + 3.417415e-01 j 
-1. 78344 le+00 + -3.417415e-01 j 


When x= 10.000000, N is 
-1.799812e-01 +-1.479897e+01 j 
8.299747e-02 +8. 161762e-02 j 


-1.799812e-01 + 1.479897e+01 j 
8.299747e-02 + -8. 161762e-02 j 


When x= 15. 000000, N is 
-7.77352 le-03 +-1.198140e-02 j 
-7.773521e-03 +1.198140e-02 j 


-2.001 004e+01 + -4.341 330e+00j 
-2.001004e+01 + 4.341330e+00 j 


When x=20.000000, N is 
-6.730543e+01 +-1.527079e+02 j 
9.609614e-03 +3.772008e-03 j 


-6.730543e+01 + 1.527079e+02 j 
9.6096 14e-03 + -3.772008e-03 j 


When x=25. 000000, N is 

-2. 1 74860e-04 +- 1 .2478 1 5e-03 j 

-2. 1 74860e-04 +1 .2478 1 5e-03 j 


-1.883595e+02 + -1.3351 19e+02j 
-1.883595e+02 + 1.33511 9e+02 j 


When x=30.000000, N is 
-1.373381e+03 +-1.286365e+03 j 
9.151 546e-04 +-2.647852e-05 j 


-1.373381e+03 + 1.286365e+03 j 
9. 1 5 1 546e-04 + 2.647852e-05 j 


When x=35. 000000, N is 
2.565397e-05 +-1.093629e-04 j 
2.565397e-05 +1.093629e-04 j 


-1.363410e+03 + -2.2 1775 9e+03 j 
-1.363410e+03 + 2.217759e+03 j 


When x=40. 000000, N is 
-1.993362e+04 +-7.270083e+03 j 
7.373917e-05 +-3.398770e-05 j 

When x=45. 000000, N is 


-1.993362e+04 + 7.270083e+03 j 
7.373917e-05 + 3.398770e-05 j 
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5.896356e-06 +-8.029830e-06 j 
5.896356e-06 +8.029830e-06 j 

When x=50. 000000, N is 
-2.389105e+05 +1.273344e+04 j 
4.833746e-06 +-5.337304e-06 j 

When x=5 5. 000000, N is 
7.602998e-07 +-4.500235e-07 j 
7.602998e-07 +4.500235e-07 j 

When x=60. 000000, N is 
-2.4218 17e+06 + 1. 188 501 e+06 j 
2.087238e-07 +-6.035369e-07 j 

When x=65. 000000, N is 
7.767776e-08 +- 1 .027 1 82e-08 j 
7.767776e-08 +1.0271 82e-08 j 

When x=70. 000000, N is 
- 1 . 9865 1 8e+07 +2 . 303671 e+07 j 
-3.961913e-09 +-5. 649675 e-08j 

When x=75. 000000, N is 
6.694335e-09 +1.863506e-09 j 
6.694335e-09 +-1.863506e-09 j 

When x=80. 000000, N is 
-1.041891e+08 +3.267910e+08 j 
-2.288082e-09 +-4.471323e-09 j 

When x=85. 000000, N is 
4.813270e-10 +3.848460e-10 j 
4.813270e-10 +-3.848460e-10j 

When x=90. 000000, N is 
3.643030e+08 +3.850378e+09 j 
-3.42171 le-10 +-2.852034e-10 j 

When x=95. 000000, N is 
2.588327e-l 1 +4.8 1 362 1 e- 1 1 j 
2.588327e-l 1 +-4.813621e-l 1 j 

When x= 100.000000, N is 
2.080449e+10 +3.832754e+10 j 
-3.783386e-l 1 +-1.136730e-ll j 



-4.335545e+03 + -2.903269e+04 j 
-4.335545e+03 + 2.903269e+04 j 



-2.389105e+05 + -1.273344e+04 j 
4.833746e-06 + 5.337304e-06 j 



8 . 3 40682e+04 + -3 . 203 1 55 e+05 j 
8.340682e+04 + 3.203 155e+05j 



-2.42181 7e+06 + - 1 . 1 8 850 1 e+06 j 
2.087238e-07 + 6.035369e-07 j 



2.281527e+06 + -2.95 3 693 e+06 j 
2.281527e+06 + 2.953693e+06 j 



-1.9865 18e+07 + -2.303671e+07j 
-3.961913e-09 + 5.649675e-08 j 



3.672609e+07 + -2.054894e+07 j 
3.672609e+07 + 2.054894e+07 j 



-1.041891e+08 + -3.267910e+08j 
-2.288082e-09 + 4.47 1 323e-09 j 



4.718089e+08 + -5.07493 1 e+07 j 
4.718089e+08 + 5.07493 1 e+07 j 



3.643030e+08 + -3.850378e+09 j 
-3.42171 le-10 + 2.852034e-10 j 



5.1 18289e+09 + 1.559857e+09 j 
5. 1 18289e+09 + -1.559857e+09 j 



2.080449e+10 + -3.832754e+10 j 
-3.783386e-l 1 + 1. 136730e-l 1 j 
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Overall, M is 



-6.842685e+04+5.333000e+05j 4.766148e+05 + -2.488562e+05 j 

4.766148e+05 +2.488562e+05 j -6.842685e+04 + -5.333000e+05 j 
Broke due to numerical inaccuracy @ E=0. 196624 eV 

The program output shows the symmetry and progression of the N matrices. At 

an interface where a step increase in potential occurs (up-step), the form of N is 

c+di a + bi 
c-di a-bi 

On the other hand, at an interface where a step decrease in potential occurs (down-step), 

N has the form 

a + bi a- bi 

N = 

1 down j- ’ 

c+ai c-di 

where, in each case, a, b, c, and d are real constants. Note that the N up and Na 0 wn matrices 
have distinctly different symmetry. It is interesting that the matrix product of all of the N 
matrices has diagonal symmetry. 

Note that, as x increases, the magnitudes of a and b, for both the N up and Nd 0W n 
matrices, get very large, while the magnitudes of c and d get very small. At the point that 
either type of N matrix’ c + di elements approach the numerical limitations of double- 
precision floating-point numbers, the determinant of the overall transfer matrix, M, 
begins to diverge from one. When this happens, the transfer matrix method breaks down, 
as it has in the above example. 
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APPENDIX B. 



TRANSFER MATRIX METHOD NUMERICAL INSTABILITY 
AS SEEN IN PROGRAM OUTPUT 



SINGLE-PRECISION FLOATING-POINT NUMBERS USED 



1. Barrier Height Varied 



% ACC withfloat 
% withfloat 
BARRIERS - 1 
BARWIDTH = 5.000000 nm 
Vo = 0.100000 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 0.200000 eV 

Broke due to numerical inaccuracy @ E=0. 010702 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 0.500000 eV 

Broke due to numerical inaccuracy @ E=0. 173308 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 1.000000 eV 

Broke due to numerical inaccuracy @ E=0. 680638 eV 



2 . 



Barrier Width Varied 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BAR WIDTH = 5.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 013768 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BAR WIDTH = 6.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.0361 12 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BAR WIDTH = 7.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 072364 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BARWIDTH = 8.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 098584 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BARWIDTH = 9.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 1 18420 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BARWIDTH = 10.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 161284 eV 
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3. 



Number of Barriers Varied 



% ACC withfloat 
% withfloat 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.013768 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 2 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 147604 eV 



% ACC withfloat 
% withfloat 
BARRIERS = 5 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.221248 eV 



B. DOUBLE-PRECISION FLOATING-POINT NUMBERS USED 



1. Barrier Height Varied 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 0.1 00000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
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Vo = 0.200000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BAR WIDTH = 5.000000 nm 
Vo = 0.500000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 1.000000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 2.000000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 5.000000 eV 

Broke due to numerical inaccuracy @ E=0. 800680 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 10.000000 eV 

Broke due to numerical inaccuracy @ E=5. 629874 eV 



2. Barrier Width Varied 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
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BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 6.000000 nm 
Vo = 0.230000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 7.000000 nm 
Vo = 0.230000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 8.000000 nm 
Vo = 0.230000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 9.000000 nm 
Vo = 0.230000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 10.000000 nm 
Vo = 0.230000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 20.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.003964 eV 
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% ACC withdouble 
% withdouble 
BARRIERS = 1 
BAR WIDTH - 30.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 107932 eV 



% ACC withdouble 
% withdouble 
BARRIERS - 1 
BARWIDTH = 40.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E— 0 .161512 eV 

% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH - 50.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 187504 eV 



3. Number of Barriers Varied 



% ACC withdouble 
% withdouble 
BARRIERS = 1 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 2 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 5 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 
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Broke due to numerical inaccuracy @ E=0. 046372 eV 



% ACC withdouble 
% withdouble 
BARRIERS - 10 
BAR WIDTH = 5.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 196624 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 15 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.2 15776 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 20 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0. 220792 eV 



% ACC withdouble 
% withdouble 
BARRIERS = 50 
BARWIDTH = 5.000000 nm 
Vo = 0.230000 eV 

Broke due to numerical inaccuracy @ E=0.227176 eV 
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APPENDIX C. 



CODE FOR C PROGRAMS WRITTEN FOR THIS THESIS 



A. TRANSFER MATRIX METHOD (methodl.c) 

/* Francis E. Spencer III */ 

/* Thesis, Summer 1997*/ 

/* Prof. Luscombe */ 

/* Inclusions */ 

#include <stdio.h> 

#include <math.h> 

#include <stdlib.h> 

/* Definitions */ 

#defme BARRIERS 1 /* dimensionless */ 

#define BARPTS 10 /* dimensionless */ 

#define BARWIDTH 5.0 /* nm */ 

#defme MAXX ((2.0*BARWIDTH*BARRIERS)+BAR WIDTH) /* nm */ 

#define XPTS ((2*BARPTS*BARRIERS)+BARPTS) /* dimensionless */ 

#define EPTS 10000 /* dimensionless */ 

#define EFFMASS 0.067 /* dimensionless */ 

#define H20VER2M 0.0381 /* eV-nm2 */ 

#define C (EFFMASS/H20VER2M) /* l/(eV-nm2) */ 

#defme Vo 0.23 /* eV */ 

#define ELEFT 0.000 1 /* eV */ 

#define ERIGHT (Vo-0.000 1) /* eV */ 

/* Function Prototypes */ 
void makeV(void); 

struct complex addc(struct complex a, struct complex b); 
struct complex subc(struct complex a, struct complex b); 
struct complex mulc(struct complex a, struct complex b); 
struct complex divc(struct complex a, struct complex b); 
struct complex expc(struct complex a); 
double absc(struct complex a); 

struct complex N1 lcomplex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N12complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex determinantcomplex(struct complex Mil, struct complex Ml 2, struct complex M21, struct 
complex M22); 

/* Global Variable Definitions */ 
struct complex 
{ 

double real; 
double imag; 
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}; 

double V[XPTS]; 
struct complex k[XPTS]; 

/* Body of Program follows: */ 

void main(void) 

/* This function controls program execution */ 

{ 

int ec,xc; 

double E,epts,xpts,T; 
struct complex 

nevvMl l,newM12,nevvM2KnewM22,oldMl l,oldM12,oldM21,oldM22,Nl l,N12,N21,N22,x,kratio,Mdet; 
x.imag=0.0; /* x is a purely real number */ 

epts=EPTS; 
xpts=XPTS; 
makeV(): 

for(ec=0.0;ec<epts;ec+=1.0) 

{ 

E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); /* DOWN-counting through E values */ 
oldMl l.imag=oldM12.imag=oldM21.imag=oldM22.imag=0.0; 
oldMl l.real=oldM22.real=1.0; 

oldM12.real=oldM21.real=0.0; /* "old 11 M initially an identity matrix */ 

for(xc=0;xc<XPTS;xc++) /* initialize k */ 

{ 

if(E>=V[xc|) /* k is purely real */ 

{ 

k[xc].real=sqrt(C*(E-V|xc])); 

k[xc].imag=0.0; 

}/* end if E>= V */ 

else /* k is purely imaginary */ 

{ 

k[xc].real=0.0; 

k [ xc ] unag= sqrt( C * ( V [ xc] -E)) ; 

}/* end else E<V */ 

}/* end for xc (init k) */ 

for(xc= 1 ;xc<XPTS;xc++) 

{ 

x . real=(xc/xpts)*M AXX; 

if((k[xc|.real!=k[xc-l].real)||(k[xc].imag!=k[xc-l).imag)) 

{ 

kratio=di vc(k [ xc] ,k [xc- 1 ] ); 

N1 1=N1 lcomplex(k[xc-l],k[xc],kratio,x); 
N12=N12complex(k[xc-l],k[xc],kratio,x); 
N21=N21complex(k[xc-l].k[xc],kratio.x); 
N22=N22complex(k[xc-l],k[xc],kratio,x); 

newMl l=addc(mulc(oldMl l,Nll),mulc(oIdM12,N21)); 
new r M12=addc(mulc(oldMl l,N12),mulc(oldM12,N22)); 
ne\vM2 1 =addc(mulc(oldM2 1 ,N 1 1 ),mulc(oldM22,N2 1 )); 
nevvM22=addc(mulc(oldM21,N12),mulc(oldM22,N22)); 

oldMl l=newMll; 
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oldM12=ne\vM12; 

oldM21=newM21; 

oldM22=nevvM22; 

}/* end if k[xc] ... */ 

}/* end for xc */ 

Mdet=determinantcomplex(newM 1 1 .newM 1 2,nevvM2 1 ,newM22); 

if((Mdet.real>l. 000 l)||(Mdet.real<0. 9999)) 

{ 

printf("Broke due to numerical inaccuracy @ E=%f eV\n\n",E); 
break; 

} /* THE "SAFETY" IS ON */ 

else 

{ 

T=sqrt(LO/(absc(nevvMl l)*absc(newMl 1))); 
printf("%. 1 2 At%. 1 2f \n",E,T); 

} 

}/* end for ec */ 

}/* end MAIN */ 

void makeV(void) 

/* This function initializes the potential energy array. V */ 

{ 

int xcounfy; 

for(xcount=0;xcount<XPTS;xcount++) 

{ 

y=xcount/BARPTS; 
if ((y % 2)==0) 

V[xcount]=0.0; 

else 

V[xcount]=Vo; 

}/* end for xcount */ 

}/* end MAKEV */ 

struct complex addc(struct complex a, struct complex b) 

/* This function adds two complex numbers passed to it */ 

{ 

struct complex sum; 
sum. real=a.real+b. real; 
sum . i mag=a . i mag+b . i mag; 
retum(sum); 

}/* end ADDC */ 

struct complex subc(struct complex a, struct complex b) 

/* This function subtracts complex number b from complex number a */ 

{ 

struct complex difference; 
difference . real=a. real -b . real ; 
difference . imag=a. imag-b. imag; 
retum(difference); 

}/* end SUBC */ 

struct complex mulc(struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 
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{ 



struct complex product; 
product.real=a.real*b.real-a.imag*b.imag; 
product.imag=a.real*b.imag+a.imag*b.real; 
retum(product); 

}/* end MULC */ 

struct complex divc(struct complex a, struct complex b) 

/* This function divides complex number a by complex number b */ 

{ 

struct complex quotient; 
double denom; 

denom=b.real*b.real+b.imag*b.imag; 
quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 
quotient. imag=(b. real * a. imag-a. real *b. imag)/deno m; 
retum(quotient); 

}/* end DIVC */ 

struct complex expc(struct complex a) 

/* This function computes the exponential of a complex number */ 

{ 

struct complex exponential; 
exponential. real=exp(a.real)*cos(a.imag); 
exponential. imag=exp(a.real)*sin(a.imag); 
retum(exponential); 

}/* end EXPC */ 

double absc(struct complex a) 

/* This function returns the magnitude of complex number a */ 

{ 

return(sqrt(a.real*a.real+a.imag*a.imag)); 

}/* end ABSC */ 

struct complex N1 lcomplex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function fmds the matrix element N1 1 of the interface */ 

{ 

struct complex one, j, N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

N=mulc(addc(one,kratio),expc(mulc(j,mulc(x,subc(kplus,kminus))))); 

N. real *=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end Nil COMPLEX */ 

struct complex N12complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function finds the matrix element N12 of the interface */ 

{ 

struct complex one, minusj, N; 
one.real=1.0; 
one.imag=0.0; 
minusj. real=0.0; 
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minusj. imag=- 1.0; 

N=mulc(subc(one,kratio),expc(mulc(minusj,mulc(x,addc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end N12COMPLEX */ 

struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function finds the matrix element N21 of the interface */ 

{ 

struct complex one, j, N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.imag=1.0, 

N=mulc(subc(one,kratio),expc(mulc(j,mulc(x,addc(kplus,kminus))))); 

N. real *=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end N21COMPLEX */ 

struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function finds the matrix element N22 of the interface */ 

{ 

struct complex one, minusj, N; 

one.real=1.0; 

one.imag=0.0; 

minusj. real=0.0, 

minusj. imag=-1.0; 

N=mulc(addc(one,kratio),expc(mulc(minusj,mulc(x,subc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

rettun(N); 

}/* end N22COMPLEX */ 

struct complex determinantcomplex(struct complex Mil, struct complex Ml 2, struct complex M21, struct 
complex M22) 

/* This function calculates the determinant of the transfer matrix as a check for accuracy */ 

{ 

retum(subc(mulc(M 1 1 ,M22),mulc(M2 1 ,M 1 2))) ; 

}/* end DETERMINANTCOMPLEX */ 



B. BACKWARD-RECURRENCE METHOD (method2.c) 

/* Francis E. Spencer III */ 

/* Thesis, Summer 1997*/ 

/* Prof. Luscombe */ 

/* Inclusions */ 

#include <stdio.h> 

#include <math.h> 
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#include <stdlib.h> 



/* Definitions */ 

#define BARRIERS 1 /* dimensionless */ 

#deline BARPTS 5000 /* dimensionless */ 

/* NOTE: set BARPTS so that there are 0.01 nm per point (BARWIDTH/BARPTS)=( 1/100) */ 
#define BAR WIDTH 5.0 /* mn */ 

#define MAXX ((2.0*BARWIDTH*BARRIERS)+BAR WIDTH) /* nm */ 

#define XPTS ((2*BARPTS*B ARRIERS)+BARPTS) /* dimensionless ♦/ 

#define EPTS 10000 /* dimensionless */ 

#define EFEMASS 0.067 /* dimensionless */ 

#define H20VER2M 0.0381 /* eV-nm2 */ 

#define C (EFFMASS/H20VER2M) /* l/(eV-nm2) */ 

#define Vo 0.23 . /* eV */ 

#define VL 0.0 /* eV */ 

#define VM 0.0 /* eV */ 

#define VR 0.0 /* eV */ 

#define ELEFT 0.000 1 /* e V */ 

#define ERIGHT (Vo-0.0001) /* eV */ 

/* Function Protot>pes */ 
void makeV(void); 

struct complex addc(struct complex a, struct complex b); 
struct complex subc(struct complex a, struct complex b); 
struct complex mulc( struct complex a, struct complex b); 
struct complex divc(struct complex a, struct complex b); 
struct complex expc(struct complex a); 
double absc(struct complex a); 

/* Global Variable Definitions */ 
struct complex 
{ 

double real; 
double imag; 

}; 

double V[XPTSJ; 

/* Body of Program follows: */ 

void main(void) 

/* This function controls program execution */ 

{ 

int ec,xc; 

double E^epts^xpts^arptsT.delta^kfkr; 

struct complex j,minusLx,m,mplusl,ro,bn[XPTS],temp,temp2; 

x.imag=0.0; /* x is a purely real number */ 

epts=EPTS; 

xpts=XPTS; 

barpts=BARPTS; 

delta=BARWIDTH/barpts; 

j.real=0.0; 

j.imag=1.0; 

minus l.real=- 1.0; 

minusl.imag=0.0; 

makeVQ; 
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for(ec=0.0;ec<epts;ec+= 1 .0) 

{ 

E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); /* DOWN-counting through E values */ 
for(xc=0;xc<XPTS;xc++) /* initialize bn*/ 

{ 

bn[xc].real=C*delta*delta*(E-V[xc])-2.0; 

bn[xc].imag=0.0; 

}/* end for xc (init k) */ 

kl=sqrt(C*(E-VL)); 
ki^sqrt(C*(E-VR)); 
mplus 1 .real=cos(kr*delta); 
mplus 1 . imag=sin(kr*delta); 

for(xc=(XPTS-BARPTS);xc>=BARPTS;xc--) /* edge of right flat zone to edge of 

left flat zone */ 

{ 

m=divc(minusl,addc(bn[xc],mplusl)); 

mplusl.real=m.real; 

mplus 1 . imag=m. iinag; 

}/* end for xc */ 

ro.real=m.real; 

ro.imag=m.imag; 

temp.real=0.0; 

temp.imag=kl*delta; 

temp2=mulc(ro,expc(temp)); 

T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(1.0- 

2.0*temp2.real+absc(ro)*absc(ro))); 

printf( ,,0 /o. 12f\t%. 12f\t%. 12f\t%. 12f\n",E,T); 

}/* end for ec */ 

}/* end MAIN */ 

void makeV(void) 

/* This function initializes the potential energy array, V (square barriers) */ 

{ 

int xcount,y; 

for(xcount=0;xcount<XPTS;xcount++) 

{ 

v=xcount/BARPTS; 

if ((y % 2)==0) 

{ 

if(xcount<BARPTS) V[xcount]=VL; /* LEFT flat zone */ 

else if(xcount>(XPTS-BARPTS)) V[xcount]=VR; /* RIGHT flat zone */ 

else V[xcount|=VM; /* flat zone between barriers */ 

}/* end if y */ 
else 

V[xcount]=Vo; /* inside a barrier */ 

}/* end for xcount */ 

}/* end MAKEV */ 

struct complex addc(struct complex a, struct complex b) 
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/* This function adds two complex numbers passed to it */ 

{ 

struct complex sum: 
sum. real=a.real+b. real: 
sum.imag=a.imag+b.imag; 
retum(sum); 

}/* end ADDC */ 

struct complex subc(struct complex a, struct complex b) 

/* This function subtracts complex number b from complex number a */ 

{ 

struct complex difference; 
difference. real=a.real-b. real; 
difference. imag=a.imag-b.imag; 
retum(difference); 

}/* end SUBC */ 

struct complex mulc(struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 

{ 

struct complex product; 
product.real=a.real*b.real-a.imag*b.imag; 
product.imag=a.real*b.imag+a.imag*b.real; 
retum(product); 

}/* end MULC */ 

struct complex divc(struct complex a, struct complex b) 

/* This function divides complex number a bv complex number b */ 

{ 

struct complex quotient; 
double denom; 

denom=b.real*b.real+b.imag*b.imag; 
quotient. real=(a.real*b.real+a.imag*b.imag)/denom; 
quotient. imag=(b.real*a.imag-a.real*b.imag)/denom; 
re tum( quotient); 

}/* end DIVC */ 

struct complex expc(struct complex a) 

/* This function computes the exponential of a complex number */ 

{ 

struct complex exponential; 
exponential. real=exp(a.real)*cos(a.imag); 
exponential. imag=exp(a.real)*sin(a.imag); 
retum(exponential); 

}/* end EXPC */ 

double absc(struct complex a) 

/* This function returns the magnitude of complex number a */ 

{ 

retum(sqrt(a.real*a.real+a.imag*a.imag)); 

}/* end ABSC */ 



C. TRANSFER MATRIX METHOD COMPARED TO ANALYTIC 
METHOD (mlwithref.c) 
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/* Francis E. Spencer III */ 

/* Thesis, Summer 1997*/ 

/* Prof. Luscombe */ 

/* Inclusions */ 

#include <stdio.h> 

#include <math.h> 

#include <stdlib.h> 

/* Definitions */ 

#define BARRIERS 1 /* dimensionless */ 

#define BARPTS 10 /* dimensionless */ 

#define BARWIDTH 5.0 /* nm */ 

#define MAXX ((2.0*BARWIDTH*BARRIERS)+BAR WIDTH) /* nm */ 

#define XPTS ((2*BARPTS*BARRIERS)+BARPTS) /* dimensionless */ 

#define EPTS 10000 /* dimensionless */ 

#define EFFMASS 0.067 /* dimensionless */ 

#define H20VER2M 0.0381 /* eV-nm2 */ 

#define C (EFFMASS/H20VER2M) /* l/(eV-nm2) */ 

#define Vo 0.23 /* eV */ 

#define ELEFT 0.000 1 /* eV */ 

#define ERIGHT (5.0* Vo) /* eV */ 

/* Function Prototypes */ 
void makeV(void); 

struct complex addc(struct complex a, struct complex b); 
struct complex subc(struct complex a, struct complex b); 
struct complex mulc(struct complex a, struct complex b); 
struct complex divc(struct complex a, struct complex b); 
struct complex expc( struct complex a); 
double absc(struct complex a); 

struct complex N1 lcomplex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N12complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex determinantcomplex(struct complex Mil, struct complex Ml 2, struct complex M21, struct 

complex M22); 

double Treffinder(void); 

/* Global Variable Definitions */ 
struct complex 
{ 

double real; 
double imag; 

}; 

double V[XPTS]; 
struct complex k[XPTS]; 

/* Body of Program follows: */ 

void main(void) 

/* This function controls program execution */ 
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{ 

int ec,xc; 

double E,epts,xpts,T; 
struct complex 

nevvMl l,newM12,ne\vM2 l,newM22,oldMl l,oldM12.oldM21,oldM22,Nl l,N12,N21,N22,x,kratio,Mdet; 
x.imag=0.0; /* x is a purely real number */ 

epts=EPTS; 
xpts=XPTS; 
makeV(); 

for(ec=O.0;ec<epts;ec+= 1.0) 

{ 

E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); /* DOWN-counting through E values */ 
oldM 1 1 . imag=oldM 1 2 . imag=oldM2 1 . imag=oldM22 . imag=0 . 0; 
oldM 1 1 . real=oldM22.real= 1 .0; 

oldM12.real=oldM21.real=0.0; /* "old" M initially an identity matrix */ 

for(xc=0;xc<XPTS;xc-H-) /* initialize k */ 

{ 

if(E>=V[xc]) /* k is purely real */ 

r 

\ 

k [ xc ] . real= sqrt ( C * (E- V [ xc] )) ; 
k[xc].imag=0.0; 

}/* end if E>=V */ 

else /* k is purely imaginary */ 

{ 

k[xc].real=0.0; 

k[xc].imag=sqrt(C*(V[xc]-E)); 

}/* end else E<V */ 

}/* end for xc (init k) */ 

for(xc= 1 ;xc<XPTS;xc++) 

{ 

x.real=(xc/xpts)*MAXX; 

if((k[xc].real!=k[xc-l].real)||(k[xc].imag!=k[xc-l].imag)) 

{ 

kratio=divc(k[xc],k[xc-l]); 

Nll=Nllcomplex(k[xc-l],k[xc],kratio,x); 

N12=N12complex(k[xc-l],k[xc],kratio v x); 

N21=N21complex(k[xc-l],k[xc],kratio.x); 

N22=N22complex(k[xc-l],k[xc],kratio,x); 

newM 1 l=addc(mulc(oldMl 1,N1 l),mulc(oldM 12,N2 1)); 
newM12=addc(mulc(oldMl l,N12),mulc(oldM12,N22)); 
newM2 l=addc(mulc(oldM2 1,N 1 l),mulc(oldM22,N2 1)); 
nevvM22=addc(mulc(oldM2 l,N12),mulc(oldM22,N22)); 

oldMll=nevvMll; 
oldM12=ne\vM12; 
oldM2 l=nevvM2 1; 
oldM22=newM22; 

}/* end if k|xc] ... */ 

}/* end for xc */ 
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/*Mdet=determinantcomplex(newM 1 1 ,newM 1 2, ne\vM2 1 , newM22); */ 

/*if((Mdet.real> 1.000 l)||(Mdet.real<0.9999)) 

{ 

printf("Broke due to numerical inaccuracy @ E=%f eV\n\n",E); 
break; 

}*/ /* THE "SAFETY 0 IS OFF */ /* restore Mdet line above, also, to turn it on again */ 

/*else 
{*/ 

T=sqrt(1.0/(absc(newMl l)*absc(newMl 1))); 

printf("%. 12At%. 12f\t%. 12f\t%. 12An H ,E,T,Tref_finder(),(T-Tref_finder())); 

/*}*/ 

}/* end for ec */ 

}/* end MAIN */ 

void makeV(void) 

/* This function initializes the potential energy array, V */ 

{ 

int xcounfy; 

for(xcount=0;xcount<XPTS;xcount-H-) 

{ 

v=xcount/BARPTS; 
if ((y % 2)==0) 

V[xcount]=0.0; 

else 

V[xcount]=Vo; 

}/* end for xcount */ 

}/* end MAKEV */ 

struct complex addc(struct complex a, struct complex b) 

/* This function adds two complex numbers passed to it */ 

{ 

struct complex sum; 
sum . real=a . real+b. real ; 
sum.imag=a.imag+b.imag; 
retum(sum); 

}/* end ADDC */ 

struct complex subc(struct complex a, struct complex b) 

/* This function subtracts complex number b from complex number a */ 

{ 

struct complex difference; 
difference. real=a. real -b. real ; 
difference . imag=a. imag-b . imag; 
retum(difiference); 

}/* end SUBC */ 

struct complex mulc( struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 

{ 

struct complex product; 

product, real =a. real *b . real -a. imag*b. imag; 

product.imag=a.real*b.imag+a.imag*b.real; 

retum(product); 
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}/* end MULC */ 



struct complex divc(struct complex a, struct complex b) 

/* This function divides complex number a by complex number b */ 

{ 

struct complex quotient; 
double denom; 

denom=b.real*b.real+b.imag*b.imag; 
quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 
quotient. imag=(b.real*a.imag-a.real*b.imag)/denom; 
retum(quotient); 

}/* end DIVC */ 

struct complex expc(struct complex a) 

/* This function computes the exponential of a complex number */ 

r 

\ 

struct complex exponential; 
exponent) al.real=exp(a.real)*cos(a.imag); 
exponential. imag=exp(a.real)*sin(a.imag); 
return(exponential); 

}/* end EXPC */ 

double absc(struct complex a) 

/* This function returns the magnitude of complex number a */ 

{ 

return(sqrt(a.real*a.real+a.imag*a.imag)); 

}/* end ABSC */ 

struct complex N1 lcomplex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function finds the matrix element N 1 1 of the interface */ 

{ 

struct complex one, j, N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

N=mulc(addc(one,kratio),expc(mulc(j,rnulc(x,subc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end Nil COMPLEX */ 

struct complex N12complex(struct complex kminus, struct complex kplus, struct complex kratio. struct 
complex x) 

/* This function finds the matrix element N12 of the interface */ 

{ 

stmet complex one, minusj. N; 

one.real=1.0; 

one.imag=0.0; 

minusj. real=0.0; 

minusj. imag=-1.0; 

N=mulc(subc(one,kratio),expc(mulc(minusj,mulc(x,addc(kplus, kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

retum(N); 



62 



}/* end N12COMPLEX */ 

struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function finds the matrix element N2 1 of the interface */ 

{ 

struct complex one, j, N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

N=mulc(subc(one,kratio),expc(mulc(j,niulc(x,addc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end N21 COMPLEX */ 

struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function finds the matrix element N22 of the interface */ 

{ 

struct complex one, minusj, N; 

one.real=1.0; 

one.imag=0.0; 

minusj. real=0.0; 

minusj. imag=-l .0; 

N=mulc(addc(one,kratio),expc(mulc(minusj,mulc(x,subc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end N22COMPLEX */ 

struct complex determinantcomplex(struct complex Ml 1, struct complex M12, struct complex M21, struct 
complex M22) 

/* This function calculates the determinant of the transfer matrix as a check for accuracy */ 

{ 

retum(subc(mulc(M 1 1 ,M22),mulc(M2 1 ,M 1 2))); 

}/* end DETERM IN ANTCOMPLEX */ 

double Tref_finder(void) 

/* This function computes the analytic transmission coefficient 
for the one-barrier system of height Vo and width BAR WIDTH */ 

{ 

struct complex num,denom,kl,k2,two,four,barwidth,j,termLterm2; 

two.real=2.0; 

two.imag=0.0; 

four.real=4.0; 

four.imag=0.0; 

barwidth.real=BAR WIDTH; 

barwi dtli . imag=0 .0 ; 

j.real=0.0; 

j.imag=1.0; 

kl.real=k[0].real; 

kl.imag=k[0].imag; 

k2.real=k[BARPTS-H].real; 

k2 . imag=k [B ARPTS+ 1 ] . imag; 
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nuxn=mulc(niulc(mulc(expc(mulc(inulc(subc(k2,k 1 ),barv\ ldth) ,j)),k2),k 1 ),four); 
term 1 =mulc(addc(k 1 ,k2 ),addc(k 1 ,k2)); 

term2=mulc(mulc(subc(kl,k2),subc(kl,k2)),expc(nmlc(mulc(mulc(k2,barvvidth)j),t\vo))); 
denom=subc( term 1 ,term2 ); 

retimi(sqrt(absc(divc(num,denom))*absc(divc(num,denom)))); 

}/* end TREF_FINDER */ 



D. BACKWARD-RECURRENCE METHOD COMPARED TO ANALYTIC 
METHOD (m2withref.c) 

/* Francis E. Spencer III */ 

/* Thesis, Summer 1997*/ 

/* Prof. Luscombe */ 

/* Inclusions */ 

#include <stdio.h> 

#include <math.h> 

#include <stdiib.h> 

/* Definitions */ 

#define BARRIERS 1 /* dimensionless */ 

#define BARPTS 500 /* dimensionless */ 

/* NOTE: set BARPTS so that there are 0.01 nm per point (BARWIDTH/BARPTS)=( 1/100) */ 

#define BAR WIDTH 5.0 /* nm */ 

#define MAXX ((2.0*BARWIDTH*BARRIERS)+BAR WIDTH) /* nm */ 

#define XPTS ((2*BARPTS*BARRIERS)+BARPTS) /* dimensionless */ 

#define EPTS 10000 /* dimensionless */ 

#define EFFMASS 0.067 /* dimensionless */ 

#define H20VER2M 0.038 1 /* eV-nm2 */ 

#define C (EFFMASS/H20VER2M) /* l/(eV-nm2) */ 

#define Vo 0.23 /* eV */ 

#define VL 0.0 /* eV */ 

#define VM 0.0 /* eV */ 

#define VR 0.0 /* eV */ 

#define ELEFT 0.0001 /* eV */ 

#define ERIGHT (5.0*Vo) /* eV */ 

/* Function Prototypes */ 
void makeV(void); 

struct complex addc(struct complex a struct complex b): 
struct complex subc(struct complex a, struct complex b); 
struct complex mulc(struct complex a, struct complex b); 
struct complex divc(struct complex a struct complex b); 
struct complex expc(struct complex a): 
double absc(struct complex a); 
double Tref_finder(double kl, double E); 

/* Global Variable Definitions */ 
struct complex 

{ 

double real; 
double imag; 

}; 

double V[XPTS]: 
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/* Body of Program follows: */ 

void main(void) 

/* This function controls program execution */ 

{ 

int ec,xc; 

double E,epts v \pts,barpts,T.delta,kl,kr; 

struct complex j,minusl,x,iumplusl,ro,bn[XPTSJ,temp,temp2; 

x.imag=0.0; /* x is a purely real number */ 

epts=EPTS; 

xpts=XPTS; 

barpts=BARPTS; 

delta=BARWIDTH/barpts; 

j.real=0.0; 

j.imag=1.0; 

minusl.real=-1.0; 

minus l.imag=0.0; 

makeV(); 

for(ec=0.0;ec<epts;ec-f-= 1.0) 

{ 

E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); /* DOWN-counting through E values */ 
for(xc=0;xc<XPTS;xc++) /* initialize bn*/ 

{ 

bn[xc].real=C*delta*delta*(E-V[xc])-2.0; 

bn[xc].imag=0.0; 

}/* end for xc (init bn) */ 

kl=sqrt(C*(E-VL)); 

ki=sqrt(C*(E-VR)); 

mplusl.real=cos(kr*delta); 

mplusl.imag=sin(kr*delta); 

for(xc=(XPTS-BARPTS);xc>=BARPTS;xc-) /* edge of right flat zone to edge of 

left flat zone */ 

{ 

m=di vc(minus 1 ,addc(bn[xc] ,mplus 1 )); 

mplusl.real=m.real; 

mplus 1 . imag=m. imag; 

}/* end for xc */ 

ro.real=m.real; 

ro.imag=m.imag; 

temp.real=0.0; 

temp . imag=kl * delta; 

temp2=mulc(ro,expc(temp)); 

T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(1.0- 

2.0*temp2.real+absc(ro)*absc(ro))); 

printf( M %.12f\t%. 12ftt% J2ftt%J2f\n^E,TTref_finder(kl,E),(T-Tref_flnder(kl,E))); 

}/* end for ec */ 

}/* end MAIN */ 

void makeV(void) 
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/* This function initializes the potential energy array, V (square barriers) */ 

{ 

int xcount,y; 

for(xcount=0;xcount<XPTS;xcount++) 

{ 

y=xcount/BARPTS; 
if ((y % 2)==0) 

{ 

if(xcount<BARPTS) V[xcount]=VL; /* LEFT flat zone */ 

else if(xcount>(XPTS-BARPTS)) V[xcount]=VR; /* RIGHT flat zone */ 

else V[xcount]=VM; /* flat zone between barriers */ 

}/* end if y */ 
else 

V[xcount]=Vo; /* inside a barrier */ 

}/* end for xcount */ 

}/* end MAKEV */ 

struct complex addc(struct complex a, struct complex b) 

/* This function adds two complex numbers passed to it */ 

{ 

struct complex suin; 
sum. real=a. real+b . real ; 
sum.imag=a.imag+b.imag; 
retum(sum); 

}/* end ADDC */ 

struct complex subc(struct complex a, struct complex b) 

/* This function subtracts complex number b from complex number a */ 

{ 

struct complex difference; 
difference.real=a.real-b.real; 
diffe re nee . imag=a . imag-b . imag; 
retum(difference); 

}/* end SUBC *1 

struct complex mulc(struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 

{ 

struct complex product; 
p roduc t . real=a . real * b . real-a. imag* b . imag; 
product. imag=a.real*b.imag+a.imag*b.real; 
retum(product); 

}/* end MULC */ 

struct complex divc(struct complex a, struct complex b) 

/* This function divides complex number a by complex number b */ 

{ 

struct complex quotient; 
double denom; 

denom=b.real*b.real+b.imag*b.imag; 

quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 

quotient.imag=(b.real*a.imag-a.real*b.imag)/denom; 

retum(quotient); 

}/* end DIVC */ 

struct complex expc(struct complex a) 
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/* This function computes the exponential of a complex number */ 

{ 

struct complex exponential; 
exponential. real=exp(a.real)*cos(a.imag); 
exponential.imag=exp(a.real)*sin(a.imag); 
retum(exponential ); 

}/* end EXPC */ 

double absc(struct complex a) 

/* This function returns the magnitude of complex number a */ 

{ 

retum(sqrt (a . real *a . real +a . i mag*a . i mag) ) ; 

}/* end ABSC */ 

double Tref_finder(double kl, double E) 

/* This function computes the analytic transmission coefficient 
for the one-barrier system of height Vo and width BAR WIDTH */ 

{ 

struct complex num,denom,kl,k2,two,four,banvidth,j,terml,term2; 

two.real=2.0; 

two.imag=0.0; 

four.real=4.0; 

four.imag=0.0; 

barwi dth . real= B AR WIDTH ; 

barwidth.imag=0.0; 

j.real=0.0; 

j.imag=1.0; 

kl.real=kl; 

kl.imag=0.0; 

if(E>=V[BARPTS+l]). /* k2 is purely real */ 

{ 

k2.real=sqrt(C*(E-V[BARPTS+l])); 

k2.imag=0.0; 

}/* end if E>= V */ 

else /* k2 is purely imaginary */ 

{ 

k2.real=0.0; 

k2.imag=sqrt(C*(V[BARPTS+l]-E)); 

}/* end else E<V */ 

num=mulc(mulc(mulc(expc(mulc(mulc(subc(k2,kl),barwidth)j)),k2Xkl),four); 
term 1 =mulc(addc(k 1 ,k2),addc(k 1 ,k2)); 

term2=mulc(mulc(subc(kl,k2),subc(kl,k2)Xexpc(mulc(mulc(mulc(k2,barwidth),j),two))); 

denom=subc(terml,term2); 

retum(sqrt(absc(divc(num,denom))*absc(divc(num,denom)))); 

}/* end TREF_FINDER */ 



E. TRANSFER MATRIX METHOD APPLIED TO PARABOLIC BARRIERS 
(mlpara.c) 

/* Francis E. Spencer III */ 

/* Thesis, Summer 1997*/ 

/* Prof Luscombe */ 

/* Inclusions */ 
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#include <stdio.h> 

#include <math.h> 

#include <stdlib.h> 

/* Definitions */ 

#deflne BARRIERS 1 /* dimensionless */ 

#define BARPTS 1000 /* dimensionless */ 

#define BARWIDTH 5.0 /* nm */ 

#defme MAXX ((2.0*BARWIDTH*BARRIERS)+BARWIDTH) /* nm */ 

#define XPTS ((2 *B ARPTS*B ARRIERS)+B ARPTS) /* dimensionless */ 

#define EPTS 1000 /* dimensionless */ 

#define EFFMASS 0.067 /* dimensionless */ 

#define H20VER2M 0.0381 /* eV-nm2 */ 

#define C (EFFMASS/H20VER2M) /* l/(eV-nm2) */ 

#define Vo 0.23 /* eV */ 

#define ELEFT 0.001 . /* eV */ 

#define ERIGHT (2.0* Vo) /* eV */ 

/* Function Prototypes */ 
void makeV(void); 

struct complex addc(struct complex a, struct complex b); 
struct complex subc(struct complex a, struct complex b); 
struct complex mulc(struct complex a, struct complex b); 
struct complex divc(struct complex a, struct complex b); 
struct complex expc(struct complex a); 
double absc(struct complex a); 

struct complex N1 lcomplex( struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N12complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x); 

struct complex determinantcomplex(struct complex Mil, struct complex Ml 2, struct complex M21, struct 
complex M2 2); 

/* Global Variable Definitions */ 
struct complex 
{ 

double real; 
double imag; 

}; 

double VfXFTSJ; 
struct complex k[XPTS), 

/* Body of Program follows: */ 

void main(void) 

/* This function controls program execution */ 

{ 

int ec,xc; 

double E,epts,xpts,T; 
struct complex 

newMl l,newM12,newM21,nevvM22,oldMl l,oldM12,oldM21,oldM22,Nl l,N12,N21,N22,x,kratio,Mdet; 
x.imag=0.(); /* x is a purely real number */ 

epts=EPTS; 
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xpts=XPTS; 

niakeVQ; 



for(ec=0;ec<EPTS;ec++) 

{ 

E=ERIGHT-(ERIGHT-ELEFT) + (ec/epts); 
oldMl l.imag=oldM12dmag=oldM21 jn^g=oldM22.imag=0.0; 
oldM 1 1 .real=oldM22 .real= 1 .0; 
oldM 1 2.real=oldM2 1 .real=0.0; 

for(xc=0;xc<XPTS;xc-H-) 

{ 

if(E>=V[xc]) 

{ 

k[xc].real=sqrt(C*(E-V[xc])); 

k[xc].imag=0.0; 

} 

else 

{ 

k[xc].real=0.0; 

k[xc].imag=sqrt(C*(V[xc]-E)); 

} 

} 

for(xc= l;xc<XPTS;xc-H+) 

{ 

x.real=(xc/xpts)*MAXX; 

if((k[xc].real!=k[xc-l].real)||(k[xc].imag!=k[xc-l].imag)) 

{ 

kratio=divc(k[xc] ,k [xc- 1 ]); 

N1 1=N1 lcomplex(k[xc-l],k[xc],kratio,x); 

N 1 2 =N 1 2complex(k [xc- 1 ] ,k [xc] ,kratio,x); 
N21=N21complex(k[xc-l],k[xc],kratio,x); 
N22=N22complex(k[xc-l],k[xc],kratio,x); 

newM 1 l=addc(mulc(oldM 1 1 ,N 1 l),mulc(oldMl 2,N2 1)); 
newM 1 2=addc(mulc(oldM 1 KN 1 2),mulc(oldM 1 2,N22»; 
newM2 l=addc(mulc(oldM2 1,N1 l),mulc(oldM22,N2 1»; 
newM22=addc(mulc(oldM2 1 ,N 12),mulc(oldM22,N22»; 

oldMll=newMll; 
oldM12=newM12; 
oldM2 1 =newM2 1 ; 
oldM22=newM22; 



} 

} 

Mdet=determinantcomplex(newMl l,newM12,newM21,newM22); 
if((Mdet.real>1.0001)||(Mdet.real<0.9999)) 

{ 

printff'Broke due to numerical inaccuracy @ E=%f eV\n\n",E); 
break; 

} 
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else 

{ 

T=sqrt(1.0/(absc(nevvMl l)*absc(newMl 1))); 
if(T>=0.001) 

printf("%.12At%. 12f\n M ,E,T); 

} 

} 

}/* end MAIN */ 
void makeV(void) 

/* This function initializes the potential energy array, V */ 

{ 

int xcountxcount2,barcount; 

double barpts,x.barV[BARPTS],tempV[XPTS],avg,pts; 

xcount=0; 

barpts=BARPTS; 

for(x=(-2.0*Vo);x<=(2.0*Vo);x+=(4.0*Vo)/barpts) 

{ 

barV[xcount]=Vo-(x*x)/(4.0*Vo); 

xcount-H-; 

}/* end for x */ 

for(xcount=0;xcount<BARPTS;xcount++) 

{ 

V[xcount]=0.0; 

}/* end for xcount */ 

for(barcount=0;barcount<BARRIERS;barcount'H-) 

{ 

xcount2=0; 

for(xcount=BARPTS+(barcount*2*BARPTS);xcount<(2*BARPTS)+(barcount*2*BARPTS);xco 

unt-H-) 

{ 

V [ xcou nt ] =bar V [ xcount2 ] ; 
xcount2++; 

}/* end for xcount */ 

for(xcount=(2*BARPTS)+(barcount*2*BARPTS);xcount<BARPTS+((barcount+l)*2*BARPTS); 

xcount-H-) 

{ 

V[xcount]=0.0; 

}/* end for xcount */ 

}/* end for barcount */ 

pts=BARPTS/ 100.0; 

for(xcoun t=0 ; xcount <XPTS ; xcount-H-) 

{ 

avg=0.0; 

if((V[xcount]!=V[\count-l])&&(xcount!=0)) 

{ 

for(xcount2=0;xcount2<(BARPTS/100);xcount2-H-) 

{ 

avg+- V[xcount+xcount2] ; 

}/* end for xcount2 (create avg) */ 

avg=avg/pts; 
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for(xcount2=0;xcount2<(BARPTS/100);xcount2-H-) 

{ 

temp V [ xcount+xcount2 ] =avg ; 

}/* end for xcount2 (write to tempV) */ 

xcount+=((B ARPTS/ 1 00)- 1 ); 

}/* end if V */ 

else /*"V+=V- M */ 

tempV[xcount]=V[xcount]; 

}/* end for xcount (create tempV) */ 

for(xcount=0;xcount<XPTS;xcount-H-) 

{ 

V[xcount]=tempV[xcount]; 

}/* end for xcount (transfer to V[xcount]) */ 

}/* end MAKEV */ 

struct complex addc(struct complex a, struct complex b) 

/* This function adds two complex numbers passed to it */ 

{ 

struct complex sum; 
sum.real=a.real+b.real; 
sum.imag=a.imag+b.imag; 
retum(sum); 

}/* end ADDC */ 

struct complex subc(struct complex a, struct complex b) 

/* This function subtracts complex number b from complex number a */ 

{ 

struct complex difference; 
difference.real=a.real-b.real; 
difference, imag^a. imag-b. imag; 
retum(difference); 

}/* end SUBC */ 

struct complex mulc(struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 

{ 

struct complex product; 

product. real=a.real*b.real-a.imag*b.imag; 

product.imag=a.real*b.imag+a.imag*b.real; 

retum(product); 

}/* endMULC */ 

struct complex divc(struct complex a, struct complex b) 

/* This function divides complex number a by complex number b */ 

{ 

struct complex quotient; 
double denom; 

denom=b.real*b.real+b.imag*b.imag; 
quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 
quoti ent. imag= (b . real *a. imag-a. real *b. i mag)/denom; 
retum(quotient); 

}/* end DIVC */ 
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struct complex expc(struct complex a) 

/* This function computes the exponential of a complex number */ 

{ 

struct complex exponential; 
exponentii.real=exp(a.real)*cos(a.imag); 
exponential. imag=exp(a.real)*sin(a.imag); 
retum(exponential); 

}/* end EXPC */ 

double absc( struct complex a) 

/* This function returns the magnitude of complex number a */ 

{ 

return(sqrt(a.real*a.real+a.imag*a.imag)); 

}/* end ABSC */ 

struct complex N1 lcomplex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function finds the matrix element N1 1 of the interface */ 

{ 

struct complex one, j, N; 

one.real=T.O; 

one.imag=0.0; 

j. real =0.0; 

j.imag=1.0; 

N=mulc(addc(one,kratio),expc(mulc(j,mulc(x,subc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end Nil COMPLEX */ 

struct complex N12complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function finds the matrix element N1 2 of the interface */ 

{ 

struct complex one, minusj, N; 

one.real=1.0; 

one.imag=0.0; 

minusj. real=0.0; 

minusj. imag=-1.0; 

N=mulc(subc(one,kratio),expc(mulc(minusj,mulc(x,addc(kplus,kminus))))); 

N.real*=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end N12COMPLEX */ 

struct complex N21complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* Tills function finds the matrix element N21 of the interface */ 

{ 

struct complex one, j, N; 

one.real=1.0; 

one.imag=0.0; 

j.real=0.0; 

j.iinag=1.0; 

N=mulc(subc(one,kratio),expc(inulc(j,inulc(x,addc(kplus,kminus))))); 
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N.real*=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end N21 COMPLEX */ 

struct complex N22complex(struct complex kminus, struct complex kplus, struct complex kratio, struct 
complex x) 

/* This function finds the matrix element N22 of the interface */ 

{ 

struct complex one, minusj, N; 

one.real=1.0; 

one.imag=0.0; 

minusj. real=0.0; 

minusj. imag=-1.0; 

N=mulc(addc(one,kratio),expc(mulc(minusj,mulc(x,subc(kplus,kminus))))); 

N. real *=0.5; 

N.imag*=0.5; 

retum(N); 

}/* end N22COMPLEX */ 

struct complex determinantcomplex(struct complex Mil, struct complex Ml 2, struct complex M21, struct 
complex M22) 

/* This function calculates the determinant of the transfer matrix as a check for accuracy */ 

{ 

retum(subc(mulc(M 1 1 ,M22),mulc(M2 1 ,M 1 2))); 

}/* end DETERMINANTCOMPLEX */ 



F. BACKWARD-RECURRENCE METHOD APPLIED TO RESONANT- 
TUNNELING DIODE (RTD) POTENTIAL (m2RTD.c) 

/* Francis E. Spencer III */ 

/* Thesis, Summer 1997*/ 

/* Prof. Luscombe */ 

/* Inclusions */ 

#include <stdio.h> 

#include <math.h> 

#include <stdlib.h> 

/* Definitions */ 

#define PI 3.141592654 
#define MAXX 50.0 
#define XPTS 5000 
#define EPTS 1000 
#define EFFMASS 0.067 
#define H20VER2M 0.0381 

#define C (EFFMASS/H20VER2M) /* l/(eV-nm2) */ 



#define Vo 1.45 


/* eV */ 


#define VL -0.1 


/* eV */ 


#define VR -0.3 


/* eV */ 


#define ELEFT (VL+0.001) 


/* eV */ 


#define ERIGHT (Vo-0.001) 


/* eV */ 



/* Function Prototypes */ 



/* dimensionless */ 

/* nm */ 

/* dimensionless */ 

/* dimensionless */ 

/* dimensionless */ 

/* eV-nm2 */ 
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void makeV(void); 

struct complex addc(struct complex a struct complex b); 
struct complex subc(struct complex a, struct complex b); 
struct complex mulc(struct complex a, struct complex b); 
struct complex divc(struct complex a struct complex b); 
struct complex expc(struct complex a); 
double absc(struct complex a); 
double absval(double x); 

/* Global Variable Definitions */ 
struct complex 
{ 

double real; 
double imag; 

}; 

double V[XPTS]; 

/* Body of Program follows: */ 

void main(void) 

/* This function controls program execution */ 

{ 

int ec,xc; 

double E,epts,xpts,barpts,T,deltakl,kr; 

struct complex j,minus 1 ,x,rrtmplus 1 ,ro,bn [XPTS],temp > temp2; 

x.imag=0.0; /* x is a purely real number */ 

epts=EPTS; 

xpts=XPTS; 

delta=MAXX/xpts; 

j.real=0.0; 

j.imag=1.0; 

minus l.real=- 1.0; 

minus l.imag=0.0; 

makeV(); 

for(ec=0.0;ec<epts;ec+= 1.0) 

{ 

E=ERIGHT-(ERIGHT-ELEFT)*(ec/epts); /* DOWN-counting through E values */ 
for(xc=0;xc<XPTS;xc++) /* initialize bn*/ 

t 

\ 

bn[xc].real == C*delta*delta*(E-V[xcl)-2.0; 

bn[xc].imag=0.0; 

}/* end for xc (init bn) */ 

kl=sqrt(C*(E-VL)); 
kr=sqrt(C*(E-VR)); 
mplusl.real=cos(kr*delta); 
mplus 1 . imag=sin(kr* delta) ; 



for(xc=(XPTS-l);xc>=0;xc— ) /* edge of right flat zone to edge of left flat zone */ 

{ 

m=divc(minus 1 ,addc(bn[xc), mplus 1 )); 
mplus l.real=m.real; 
mplus 1 . imag=m. imag; 

}/* end for xc */ 

ro.real=m.real; 
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ro.imag=m.imag; 

temp.real=0.0; 

teinp.imag=kl*delta; 

tenip2=mulc(ro,expc(temp)); 

T=sqrt((sin(kl*delta)/sin(kr*delta))*(4.0*sin(kl*delta)*ro.imag)/(1.0- 

2.0*temp2.real+absc(ro)*absc(ro))); 

printf( M %f\t%f\n , ',E,T); 

}/* end for ec */ 

}/* end MAIN */ 

void makeV(void) 

/* This function initializes the potential energy array, V */ 

{ 

int xcounf xcount2,PTS; 

double barpts,x,xt\vo,tempV[XPTS),avg,pts,xpts,last; 
xpts=XPTS; 

for(xcount=0;xcount<XPTS;xcount*H-) 

{ 

x=xcount*(PI/xpts)-PI/2.0; 

V[xcount]=(-0. l*atan(x))-0.2; 

}/* end for xeount (create the AT AN shape) */ 

for(xcount=0;xcount<XPTS;xcount-H-) 

{ 

x=xcount*(MAXX/xpts); 
xtwo=xcount*(PI/xpts)-PI/2.0; 
if((x>= 1 8. 8)&&(x< 19.0)) 

V[xcount|=((-0.1*atan(xtwo))-0.2)+(1.55/0.2)*(x-18.8); 
else if((x>=19.0)&&(x<21.0)) 

V[xcount]= 1 .45-(0.05/2.0)*(x- 1 9.0); 
else if((x>=2 1 .0)&&(x<2 1 . 2)) 

V[xcount]=l .40-( 1 .55/0.2)*(x-2 1 .0); 
else if((x>=21.8)&&(x<22.0)) 

V[xcount]=((-0. 1 *atan(xtwo))-0.2)-(0. 2/0.2)*(x-2 1.8); 
else if((x>=22.0)&&(x<24.0)) 

V[xcount]=-0.35-(0.05/2.0)*(x-22.0); 
else if((x>=24.0)&&(x<24.2)) 

V[xcount]=-0.40+(0.2/0.2)*(x-24.0); 
else if((x>=24.8)&&(x<25.0)) 

V[xcount]=((-0. 1 *atan(xtwo))-0.2)+( 1 .55/0. 2)*(x-24 .8); 
else if((x>=25.0)&&(x<27.0)) 

V[xcountJ=1.35-(0.05/2.0)*(x-25.0); 
else if((x>=27.0)&&(x<27.2)) 

V[xcount|= 1 .30-( 1 .55/0.2)*(x-27.0); 
else; /* leave V unchanged */ 

}/* end for xeount (add the 3 peaks) */ 

}/* end MAKEV */ 

double absval(double x) 

/* This function returns the absolute value of a double, x. */ 

{ 
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if(x<0) x=-x; 
retum(x); 

}/* end ABSVAL */ 

struct complex addc(struct complex a, struct complex b) 

/* This function adds two complex numbers passed to it */ 

{ 

struct complex sum; 
sum.real=a.real+b.real; 
sum.imag=a.imag+b. imag; 
retum(sum); 

}/* end ADDC */ 

struct complex subc(struct complex a, struct complex b) 

/* This function subtracts complex number b from complex number a */ 

{ 

struct complex difference; 
difference, real =a. real -b . real ; 
difference. imag=a. imag-b . imag; 
retum(difference) ; 

}/* end SUBC */ 

struct complex mulc(struct complex a, struct complex b) 

/* This function multiplies two complex numbers passed to it */ 

{ 

struct complex product; 
product.real=a.real*b.real-a.imag*b.imag; 
product.imag=a.real*b.imag+a.imag*b.real; 
retum(product); 

}/* end MULC */ 

struct complex divc(struct complex a, struct complex b) 

/* This function divides complex number a bv complex number b */ 

{ 

struct complex quotient; 
double denom; 

denom=b.real*b.real+b.imag*b.imag; 
quotient.real=(a.real*b.real+a.imag*b.imag)/denom; 
quotient. imag=(b.real*a.imag-a.real*b.imag)/denom; 
retum(quotient); 

}/* end DIVC */ 

struct complex expc(struct complex a) 

/* This function computes the exponential of a complex number */ 

{ 

struct complex exponential; 

expo nential . real=exp(a. real) * cos(a. imag) ; 

exponential . imag=exp(a. real) *sin(a. imag); 

retum(exponential); 

}/* end EXPC */ 

double absc(struct complex a) 

/* This function returns the magnitude of complex number a */ 

{ 

retum(sqrt(a.real*a.reaH-a.imag*a.iniag)); }/* end ABSC */ 
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APPENDIX D 



FIGURES 



1 Square Barrier, 0.23 eV, 5.0 nm 
Method 1 (Transfer Matrices) 




E (eV) 



Figure 1. Transfer Matrix Method Applied to a Single Square Barrier 
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1 Square Barrier, 0.23 eV, 5.0 nm 

Method 1 (Transfer Matrices, "Safety" Off) 

minus Analytic Expression (for 1 square barrier) 
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Figure 2. Difference Between Analytic Solution and Transfer Matrix Solution, for the 

Case of One Square Barrier 
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Three Parabolic Barriers 
(Partitioned Into Sub-Bins) 




Point Index (dimensionless) 



Figure 3. Potential Energy Profile for Three Parabolic Barriers 
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2 Parabolic Barriers, 0,23 eV, 5,0 nm 
Method 1 (Transfer Matrices) 




E (eV) 



Figure 4. Transfer Matrix Method Applied to a Two-Parabolic-Barrier Potential 
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3 Parabolic Barriers, 0,23 eV, 5 . 0 nm 
Method 1 (Transfer Matrices) 




Figure 5. Transfer Matrix Method Applied to a Three-Parabolic-Barrier Potential 
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(dimensionless) 



& Parabolic Barriers , 0.23 eV, 5.0 nm 
Method 1 (Transfer Matrices) 




E (eV) 



Figure 6. Transfer Matrix Method Applied to a Five-Parabolic-Barrier Potential 
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Potential Energy versus Distance for an 
AlAs/InGaAs/InAs Resonant Tunneling Diode, 
as in Luscombe's Article 




Figure 7. Potential Energy Profile for a Resonant Tunneling Diode 
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RTD Potential Profile (erroneous) 
TRANSFER MATRIX " SAFETY ■ OFF 




E (eV) 



Figure 8 Transfer Matrix Method Applied to the Resonant Tunneling Diode Potential 
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1 Square Barrier, 0.23 eV, 5.0 nm 
Method 2 (Backward Recurrence) 
Compared to Analytic Expression 




E (eV) 



Figure 9. Backward-Recurrence Method Applied to the Square Barrier of Figure 1 
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1 Square Barrier, 0,23 eV, 5.0 ran 
Method 2 (Backward Recurrence) 
Minus Analytic Expression 




E (eV) 



Figure 10. Difference Between Analytic Solution and Backward-Recurrence Solution, 

for the Case of One Square Barrier 
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2 Parabolic Barriers, 0.23 eV, S.O nm 
Method 2 (Backward Recurrence) 




Figure 1 1 . Backward-Recurrence Method Applied to the Two-Parabolic-Barrier 

Potential of Figure 4 
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3 Parabolic Barriers, 0.23 eV, 5.0 nm 
Method 2 (Backward Recurrence) 




Figure 12. Backward-Recurrence Method Applied to a Three-Parabolic-Barrier Potential 
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5 Parabolic Barriers, 0.23 eV, 5.0 nm 
Method 2 (Backward Recurrence) 




E (eV) 



Figure 13. Backward-Recurrence Method Applied to a Five-Parabolic-Barrier Potential 



5 Parabolic Barriers, 0.23 eV, 5.0 run 
Method 2 (Backward Recurrence) 

Using 10,000 Energy Points 




Figure 14 Backward-Recurrence Method Applied to the Five-Parabolic-Barrier Potential 
of Figure 13, Using Ten Thousand Energy Points 
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RTD Potential Profile 

Using Method 2 (backward recurrence) 

(10,000 E points ) 




E (eV) 



Figure 15 Backward-Recurrence Method Applied to the Resonant Tunneling Diode 

Potential 
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